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ABSTRACT 


This  dissertation  provides  detailed  analysis,  performance  evaluation,  and 
modeling  of  a  source  localization  method  that  is  based  on  medium  calibration. 

In  the  past,  processing  the  signals  received  on  an  array  of  sensors  in  the  ocean  to 
locate  an  acoustic  source  has  been  treated  by  methods  which  are  based  on 
detailed  knowledge  of  the  sound  speed  structure  between  the  source  and  the 
receiving  array. 

The  general  problem  considers  an  array  of  sensors  distributed  throughout 
the  water  column.  From  the  signal  measured  at  this  array,  the  location  (range 
and  depth)  of  the  unknown  source  is  determined.  A  holographic  array  processing 
(HAP)  algorithm  is  developed,  which  uses  measured  field  data  from  the 
unknown  source,  and  measurements  from  a  reference  source  deployed  near  the 
target  and  transmitting  at  several  depths.  A  virtual  image  of  the  wavefront  from 
the  unknown  source  is  constructed  at  the  location  of  the  reference  array  by  pro¬ 
cessing  the  measurements  at  the  receiving  array. 

Conventional  array  processing  algorithms,  such  as  matched  field  processing 
(MFP),  require  precise  knowledge  of  the  medium  between  the  source  and  the 
receiving  array,  but  the  HAP  relaxes  this  stiff  requirement.  It  calibrates  the  inte¬ 
grated  effect  of  a  great  portion  of  the  medium,  and  geoacoustic  parameter  esti¬ 
mation  is  needed  only  for  a  small  portion  of  the  ocean,  between  the  unknown 
source  and  the  virtual  array.  The  virtual  array  is  constructed  by  moving  a  refer¬ 
ence  source  to  incremental  depths  of  the  water  column. 

The  detection  criterion  being  used  is  based  on  the  maximum  intensity  of  the 
reconstructed  field,  which  can  be  generated  by  weighting  a  replica  field  from  a 
test  source  at  a  hypothetical  location  with  the  filtered  wavefront  from  the 
unknown  source  at  the  virtual  array. 

Theoretical  analysis  is  provided  using  the  WKB  approximation  for  a  range 
dependent  ocean  with  both  N2-bilinear  and  exponential  stratification  of  the 
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sound  speed  profile  (SSP).  Also,  numerical  simulation  is  performed  using  a 
wide  angle  PE  code1  (FEPE)  for  a  range  dependent  analytical  SSP,  and  real 
sound  speed  data  from  the  North  Pacific  ocean.  The  results  of  the  analysis  and 
simulation  show  the  possibility  of  localizing  a  source  at  large  distance  with 

great  accuracy. 

In  deep  water,  the  ideal  situation  of  having  a  receiving  array  and  reference 
source  that  span  the  water  column  is  prohibitively  complex.  How  truncation  of 
both  the  receiving  and  virtual  arrays  affects  the  performance  of  the  HAP  algo¬ 
rithm  is  investigated.  Simulation  for  a  range  dependent,  horizontally  stratified 
ocean  is  performed  using  FEPE.  It  is  shown  that  the  holographic  array  process¬ 
ing  algorithm  is  more  susceptible  to  the  length  of  the  virtual  array  than  the 
length  of  the  receiving  array. 

Detailed  analysis  is  provided  for  the  case  where  the  signals  at  the  receiving 
array  are  observed  in  the  presence  of  additive,  spatially  correlated  ambient 
noise.  For  the  numerical  simulation,  a  simplified  but  representative,  isotropic 
ambient  noise  model  is  used.  It  is  shown  that  the  power  level  of  both  the 
unknown  and  the  reference  source  is  critical  to  the  performance  of  the  HAP 
algorithm  in  the  presence  of  ambient  noise. 
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1  PE  code  isa  numerical  model  for  the  wave  propagation  in  a  range  dependent 
medium  based  on  the  parabolic  approximation  to  the  wave  equation  solution. 
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INTRODUCTION 


Processing  the  signals  received  on  an  array  of  sensors  in  the  ocean  to  locate 
a  sound  source  has  been  treated  by  methods  which  are  based  on  detailed  knowl¬ 
edge  of  the  sound  speed  structure  between  the  source  and  the  receiving  array. 

This  dissertation  provides  detailed  analysis,  performance  evaluation,  and 
modeling  of  a  source  localization  method  that  is  based  on  medium  calibration. 

1.0  Background 

Underwater  acoustics  comprises  a  major  part  of  ocean  exploration  research. 
Only  sound  waves  can  propagate  in  water  over  large  distances.  Practically  all 
kinds  of  telemetry,  communication,  localization,  and  remote  sensing  of  water 
masses  and  the  ocean  bottom  use  sound  waves.  Propagating  over  thousands  of 
kilometers  in  the  ocean,  they  bring  information  on  remote  objects,  earthquakes, 
eruptions  of  volcanoes,  distant  storms,  and  global  warming. 

The  most  characteristic  feature  of  the  ocean  medium  is  its  inhomogeneous 
nature.  There  are  two  kinds  of  inhomogeneities,  regular  and  random,  and  both 
strongly  influence  the  underwater  sound  field.  For  example  the  regular  variation 
of  the  sound  velocity  with  depth  leads  to  the  formation  of  the  underwater  sound 
channel  (USC),  and  as  a  consequence,  to  long  range  propagation.  The  random 
inhomogeneities  give  rise  to  scattering  of  sound  waves  and  therefore  to  fluctua¬ 
tions  in  the  sound  field. 

The  general  problem  we  consider  here  involves  an  array  of  sensors  distrib¬ 
uted  throughout  the  water  column.  From  the  signal  measured  at  this  array,  the 
location  (range  and  depth)  of  the  unknown  source  is  determined.  A  holographic 
array  processing  (HAP)  algorithm  is  developed,  which  uses  measured  field  data 
from  the  unknown  source  and  measurements  from  a  reference  source  deployed 
near  the  target  and  transmitting  at  several  depths.  A  virtual  image  of  the  wave- 
front  from  the  unknown  source  is  constructed  at  the  location  of  the  reference 
array  by  processing  the  measurements  at  the  receiving  array. 

Conventional  array  processing  algorithms,  such  as  matched  field  processing 
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(MFP),  require  precise  knowledge  of  the  medium  between  the  source  and  the 
receiving  array,  but  the  HAP  relaxes  this  stiff  requirement.  It  calibrates  the  inte¬ 
grated  effect  of  a  great  portion  of  the  medium,  and  geoacoustic  parameter  esti¬ 
mation  is  needed  only  for  a  small  portion  of  the  ocean,  between  the  unknown 
source  and  the  virtual  array. 

The  detection  criterion  being  used  is  based  on  the  maximum  intensity  of  the 
reconstructed  field,  which  can  be  generated  by  weighting  a  replica  field  from  a 
test  source  at  a  hypothetical  location,  with  the  filtered  wavefront  from  the 
unknown  source  at  the  virtual  array.  It  will  be  demonstrated  that  the  HAP  algo¬ 
rithm  satisfies  the  detection  criterion  when  the  test  source  is  at  the  location  of 
the  unknown  target. 

Our  prime  concern  in  this  thesis  will  be  to  pose  the  HAP  algorithm  in  pre¬ 
cise  terms  and  evaluate  its  performance  for  various  system  parameters  and  oce¬ 
anic  environments. 

2.0  Approach 

The  analytical  tools  used  in  this  research  come  from  the  fields  of  normal 
mode  theory,  multidimensional  geometry,  and  multivariate  statistics.  The  ocean 
model  is  assumed  to  be  cylindrically  symmetric  with  pressure  release  surface 
and  absorbing  bottom.  Advantage  is  taken  of  the  WKB  approximation  to  the 
wave  equation  solution  in  a  range  dependent,  horizontally  stratified  medium. 

The  computer  simulation  tools  used  in  this  thesis  are  based  on  a  wide  angle 
PE  code,  FEPE1,  and  a  normal  mode  code,  for  an  adiabatic  medium,  based  on 
the  WKB  approximation2.  The  geoacoustic  data  being  used  is  both  analytical 
and  experimental. 

The  first  chapter  will  shed  some  light  on  the  intricacy  of  underwater  wave 
propagation.  The  structure  of  the  energy  flow  in  the  ocean  makes  the  source 
localization  problem  very  complex.  Readers  who  are  familiar  with  the  funda- 


1.  FEPE  was  written  by  Mike  Collins.  The  author  modified  the  code  for  the  HAP  application. 

2.  This  code  was  developed  by  the  author  using  Matlab  software. 
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mentals  of  underwater  wave  propagation  can  proceed  to  chapter  2. 

Chapter  2  describes  the  problems  associated  with  conventional  source  local¬ 
ization  techniques  and  presents  the  motivation  behind  the  HAP  method.  The 
HAP  method  has  its  genesis  in  optics,  where  the  need  to  improve  the  resolution 
of  imaging  systems  had  led  Goodman,  Huntley,  Jackson,  and  Lehman  [1]  to  pro¬ 
pose  the  use  of  a  reference  source  in  the  vicinity  of  the  object,  record  on  a  holo¬ 
gram  the  intensity  caused  by  both  the  object  and  the  reference,  and  then  expose 
the  hologram  to  a  coherent  reference  only.  This  results  in  an  improved  image  of 
the  object.  This  chapter  presents  the  HAP  algorithm  for  source  localization  in  a 
range  dependent  medium. 

Chapter  3  provides  theoretical  and  numerical  evaluation  of  the  performance 
of  the  HAP  algorithm.  The  theoretical  analysis  is  provided  for  an  arbitrary  adia¬ 
batic  medium  by  the  WKB  approximation  to  the  wave  equation.  Detailed  results 

r\ 

are  obtained  for  a  range  dependent  ocean  with  both  N  -bilinear  and  expo¬ 
nential  stratification  of  the  sound  speed  profile  (SSP).  Also,  numerical 
simulation  is  performed  using  a  wide  angle  PE  code,  FEPE,  for  a  range 
dependent  analytical  SSP,  and  real  sound  speed  data  from  the  North 
Pacific  ocean. 

Chapter  4  presents  the  effect  of  truncating  the  receiving  and  virtual  array  on 
the  performance  of  the  HAP  algorithm.  It  will  be  demonstrated  that  the  perfor¬ 
mance  of  the  HAP  algorithm  is  more  severely  degraded  by  truncating  the  virtual 
array  than  by  truncating  the  receiving  array. 

Chapter  5  provides  detailed  analysis  for  the  performance  of  the  HAP  algo¬ 
rithm  in  the  presence  of  additive,  spatially  correlated  ambient  noise. 
Numerical  simulation  is  provided  for  the  case  when  both  the  target  and  refer¬ 
ence  signals  are  contaminated  by  the  ambient  noise.  A  simplified,  but  represen¬ 
tative  isotropic  ambient  noise  model  will  be  used  for  the  numerical  simulation. 

3.0  Contribution 


There  are  four  areas  within  which  contributions  have  been  made: 
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First,  the  geometry  of  the  emitter  location  problem  is  defined.  The  distribu¬ 
tion  of  the  receiving  and  reference  array  relative  to  the  location  of  the  unknown 
source  is  determined.  The  issue  of  sampling  the  receiving  and  virtual  array  is 
discussed. 

Second,  the  holographic  array  processing  method  is  implemented  for  source 
localization  in  a  range  dependent  ocean.  Detailed  theoretical  and  numerical 
analysis  is  provided.  The  theoretical  analysis  is  based  on  the  WKB  approxima¬ 
tion  to  the  solution  of  the  wave  equation.  The  numerical  modeling  is  performed 
using  normal  mode  code  for  an  adiabatic  medium  and  a  high-order  PE  code 
(FEPE)  which  is  based  on  the  Pade  series  approximation  to  the  wave  equation. 
The  performance  of  the  holographic  method  is  compared  with  the  performance 
of  conventional  methods,  such  as  matched  field  processing. 

Third,  the  performance  of  the  holographic  method  in  different  environmen¬ 
tal  conditions  is  evaluated.  Detailed  theoretical  and  numerical  analysis  is  pro¬ 
vided  for  the  effect  of  the  sound  speed  profile  variation  with  range  and  the 
presence  of  ambient  noise. 

Fourth,  the  performance  of  the  holographic  method  with  varying  system 
parameters  is  evaluated.  Detailed  theoretical  and  numerical  analysis  is  provided 
for  the  effect  of  implementing  a  truncated  receiving  and  a  truncated  virtual 
array. 
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4.0  Glossary 

a  =  (r',  z')  The  test  source  location 
C  (r)  Eq.  (3-11) 

c(z,r)  Sound  speed  [m/sec]  at  depth  z  and  range  r 

c,0  Axial  sound  speed  at  range  zero  (reference  sound  speed) 

c4  Axial  sound  speed  [m/sec] 

clp  Sound  speed  at  a  turning  point 

ETR  (r',  z')  The  error  term  due  to  truncating  the  receiving  array 

E1^  (r',  z ')  The  error  term  due  to  truncating  the  virtual  array 

FEPE  Computer  code  based  on  higher-order  approximation  to  the  PE 

f  The  acoustical  frequency  in  [Hz] 

G(r,z)  Green’s  function 

Gn(d,f)  Spatial  cross-spectral  density 

HAP  Holographic  Array  Processing 

H(zr)  The  holographic  image  at  reference  depth  zR 

h  The  water  depth  at  the  receiving  array 

hR  The  water  depth  at  the  reference  array 

h  The  length  of  the  truncated  receiving  array 

hR  The  length  of  the  truncated  virtual  array 


I 


Intensity 


6 


k 

MFP 

N(z,r) 

N(  6) 

PE 

P(r,z) 

P, 

Pr 

Q(r',z') 
^12  (d,  t) 
S  or  s 

S„ 

Slo 

SNR 

SSP 

TL 

z 

zitp 

zutp 

a 

a  (m) 


2n  co 
wavenumber  =  -*-  =  — 

A.  c 

Matched  Field  Processing 
Normalized  index  of  refraction 
Ambient  noise  intensity  per  solid  angle 
Parabolic  Equation 

Complex  pressure  field  at  range  r  and  depth  z 
Pressure  field  from  the  target 
Pressure  field  from  the  reference  source 
The  reconstructed  field  at  range  r'  and  depth  z' 
Temporal-spatial  cross  correlation  function 
The  horizontal  gradient  in  the  sound  speed 
The  slope  of  the  upper  portion  of  the  N2-bilinear  profile 
The  slope  of  the  lower  portion  of  the  N2-bilinear  profile 
Signal  to  noise  ratio 
Sound  Speed  Profile 
Transmission  Loss 
Depth  [m] 

Depth  of  the  lower  turning  point 
Depth  of  the  upper  turning  point 
Attenuation  factor  [dB/m] 

Eq.  (3-15) 
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P(r)  Eq.  (3-16) 

7  The  vertical  eigennumber 

A  =  r,~rR  The  separation  between  the  target  and  the  reference  array 
A  =  r'  -rR  The  separation  between  the  test  source  and  the  reference  array 
S£  =  r'  —  rt  The  separation  between  the  test  source  and  the  target 
Eq.  (3-1),  the  phase  term  in  the  modal  expression  of  Q 
0(z1,z2)  The  phase  accumulation  in  the  vertical  between  zx  and  z2 
X  The  acoustic  wave  length 

£  The  horizontal  eigennumber  (eigenvalue) 

4m  The  modal  eigennumber 

a2  The  power  of  the  unknown  target  [dB] 

a2  The  power  of  the  reference  source  [dB] 

c2t  The  intensity  of  the  ambient  noise  [dB] 

4>m  The  exit  angle  that  corresponds  to  mode  m 

^(r.z)  The  eigenfunction  at  range  r  and  depth  z 

^  The  modal  eigenfunction 

Q  Solid  angle 

co  The  acoustical  frequency  in  radians;  co  =  2nf 


CHAPTER  1 


Underwater  Wave  Propagation 

The  ocean,  together  with  its  boundaries,  forms  a  remarkably  complex 
medium  for  the  propagation  of  sound.  Propagating  through  the  ocean,  an  under¬ 
water  sound  signal  becomes  delayed,  distorted,  smeared,  and  weakened. 

1-1  Transmission  Loss 

Transmission  loss  quantitatively  describes  the  weakening  of  signal  between 
a  reference  point  (one  meter  from  the  source)  and  the  point  under  consideration. 
If  I0  is  the  intensity  at  the  reference  point  and  Ij  is  the  intensity  at  a  distant 
point,  then  transmission  loss  at  the  distant  point  is 


TL  =  lOlog 


(1-1) 


Transmission  loss  summarizes  the  effect  of  a  variety  of  underwater  propaga¬ 
tion  phenomena,  and  can  be  considered  as  the  sum  of  a  loss  due  to  spreading  and 
a  loss  due  to  attenuation. 

1-1.1  Spreading  Loss 

Spreading  loss  is  a  geometrical  effect  representing  the  weakening  of  the  sig¬ 
nal  as  it  spreads  outward  from  the  source.  Two  kinds  of  spreading  are  of  special 
interest  in  underwater  propagation. 

•  Spherical  spreading:  The  power  generated  by  the  source  is  radiated  equally  in 
all  directions,  and,  assuming  no  loss  in  the  medium,  then  the  power  crossing 
all  spheres  must  be  the  same: 


P  =  4?t  rll0  =  4tc  rJ/j  =  ...  =  4ur2nIn 


(1-2) 


where  tq  is  the  reference  distance  (usually  one  meter),  and  In  is  the  intensity  of 
sound  at  range  rn.  Transmission  loss  in  decibels  is 
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TL  =  lOlog^j  =  lOlog^jp  =  201og  (rn)  .  (1-3) 

•  Cylindrical  spreading:  The  power  radiated  by  the  source  is  distributed  over 
the  surface  of  a  cylinder  having  a  radius  equal  to  the  range  and  a  height  equal 
to  the  distance  between  the  upper  and  lower  boundaries  of  the  medium: 

P  =  2nr0HI0  =  2nrxHIx  =  ...  =  2jc  rnHIn  (14) 

and  the  transmission  loss  in  decibels  is 


FIGURE  1-1.  (a)  Spherical  spreading  in  unbounded  medium,  (b)  Cylindrical  spreading 
between  two  parallel  planes. 
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1-1.2  Attenuation  Loss 

Attenuation  loss  includes  the  effects  of  absorption,  scattering,  and  leakage 
out  of  sound  channels.  The  presence  of  attenuation  will  lead  to  a  modification  to 
the  transmission  loss  expressions.  For  spherical  spreading 

TL  =  201og  (r)  +  ar  ,  (W 

and  for  cylindrical  spreading 


TL  =  lOlog  (r)  +  ar  ,  '•w; 

where  r  is  the  range  in  meters,  and  a  is  the  attenuation  in  dB/tneter. 

Extensive  theoretical  and  experimental  analysis  have  been  conducted  [21  to 
study  the  behavior  and  causes  of  wave  attenuation  in  the  sea  water.  It  has  been 
found  that  the  main  factors  that  cause  the  absorption  of  sound  in  the  sea  water 

are: 

.  volume  viscosity  due  to  some  relaxation  processes  in  sea  water, 

.  scattering  of  sound  waves  by  inhomogeneities  of  different  kind. 

The  empirical  formula,  due  to  Marsh  and  Schulkin  [31,  for  attenuation  coef 
ficient  in  sea  water  at  frequencies  between  3  kHz  and  0.5  MHz  ts 


a  =  8.68  x  10J 


x2.34  x  10  6S//  3.38  x  10  6/ 


/r 


(1-6.54x10 ~V) 


(1-8) 


— *  - - — > 

1<P<1000,  f  is  the  frequency  [kHz],  fT  -  21.*  x  iu 

a, ion  frequency  [kHz],  and  T  is  temperature  [°C],  Since  T  ranges  between  0  and 
30°,  then  fj  ranges  between  59  and  210  kHz. 

At  low  frequencies  <100  Hz  •  3  kHz)  the  attenuation  is  better  described  by 
Thorp’s  formula  [4] 
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Q.ll/  +  44/  >  (1-9) 

1+/  4100+/ 


where  a  is  the  attenuation  coefficient  in  dB/km. 


»re  a  is  tne  aucnuauu..  - 

a  as  a  function  of  frequency  on 
c.  1  9  chows  the  attenuation  coefficient ,  a,  a 
Fig.  1-2  snows  ,  f  T— 20°C  S=35  parts  per  thou- 

a  logarithmic  scale.  The  curve  is  calculated  for  T-20  C, 
sand,  and  P=100  kg/cm2. 


1-2  Sound  Speed  Structure  in  the  Ocean 

,  fSSP1  reveals  the  variation  of  sound  velocity  wit 

The  sound  speed  profile  (SSP)  function  of  depth, 

depth.  Several  empirical  expressions  for  the  soun  accurate 

temperature,  and  salinity  were  proposed  by  many  researchers. 
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formula  is  the  one  developed  by  Wilson  [5],  However,  it  is  a  complicated 
expression  containing  many  high-order  and  cross-product 

•  ,  without  the  use  of  a  computer.  Leroy  [6]  proposed  simpler  expires 

‘sTons  that  are  believed  to  fit  the  data  as  well  as  the  Wilson  formula.  Leroy  s  or- 
mula,  often  sufficient  for  practical  work,  is 


c=  1492.9 +  3  (T- 10) -0.006  (T- 10) 2 -0.04  (T  18)  + 

Z 

1.2  (S  -  35)  -  0.01  (T -  18)  (S  -  35)  +  ^ 


(1-10) 


where  c  is  the  sound  speed,  T  is  temperature  [°C],  S  is  the  salinity  [parts  per 
thousand]  (0<S<37),  and  Z  is  the  depth  in  meters. 

1-2.1  Layers  of  the  SSP  in  the  Ocean 

A  tvnical  deep  water  SSP  is  shown  in  Fig.  1-3.  The  profile  may  be  divided 
al  lavers  having  different  thicknesses,  characteristics  and  occurrence 
SZTJT. UL  season,  time  of  the  day,  and  meteorological  condi- 

tions.  These  layers  are: 

.  Surface  layer,  just  below  the  sea  surface.  The  sound  speed  in  this  layer  ts  sus¬ 
ceptible  to  daily  and  local  changes  of  heating,  cooling  and  wind  action 

layer  may  contain  a  mixedT^  £" “ -d^o  bJtrapped 
action  of  Under  prolonged  calm  and  sunny  conditions 

lUxed  layer  disappears  and  is  replaced  by  water  in  which  the  temperature 
decreases  with  depth. 

nr : = sssssss  '=•  - 
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Main  therm ocline,  just  below  the  seasonal  thermocline.  This  layer  is  affected 
only  slightly  by  seasonal  changes.  Temperature  and  velocity  decreases  with 

depth. 

Isothermal  layer,  below  the  main  thermocline.  This  layer  has  nearly  constant 
temperature  near  4°C.  The  speed  of  sound  increases  with  depth  because  of 
effect  of  pressure  (depth)  on  the  sound  speed. 

SOFAR  (Sound  Fixing  And  Ranging)  axis.  Between  the  negative  sound  speed 
gradient  of  the  main  thermocline  and  the  positive  sound  speed  gradient  of  the 
isothermal  layer,  there  is  a  minimum  sound  speed.  Sound  traveling  a.  great 
depths  tends  to  be  bent  or  focused  by  refraction  toward  the  SOFAR  axis. 


sound  speed  (m/sec) 
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1-3  Patterns  of  Sound  Propagation  in  the  Ocean 

As  the  sound  wave  leaves  the  source,  it  follows  different  ray  paths,  and  the 
wave  propagates  along  each  path  with  different  velocity,  making  the  wavefront 
nonplanar.  Two  major  types  of  paths  exist,  reflected  and  refracted  paths. 

Reflected  paths  are  created  upon  reflection  from  the  surface  and/or  the  bot¬ 
tom;  these  paths  are  not  important  for  long  range  propagation,  because  of  the 
attenuation  (loss)  associated  with  them  upon  reflection.  Refracted  paths,  on  the 
other  hand,  carry  with  them  the  relative  information  about  the  source  and  vol¬ 
ume  of  the  medium.  These  paths  are  confined  between  turning  points  that  can  be 

found  from  Snell’s  law 


cos(pa  costp,.  coscpN 

~a  ci  CN 

where  <pa  is  the  initial  ray  exit  angle,  and 


(1-11) 


i*n<f 

(Pi  is  the  ray  angle  with  the  horizontal  at  depth  zj, 

(Pjp  is  the  array  angle  with  the  horizontal  at  the  turning  point,  <Ptp— 0  > 

cj  is  the  sound  speed  at  depth  Zj, 

ctp  is  the  sound  speed  at  the  turning  point. 

The  sound  energy  is  trapped  in  the  USC  (Underwater  Sound  Channel)  and 
does  not  reach  the  boundaries;  therefore,  the  wave  will  not  undergo  scattering 
and  absorption  at  these  boundaries.  The  greater  Ac  =  ctp-ca  ,  the  larger  is  the 
interval  of  angles  in  which  the  rays  are  trapped;  i.e.  the  waveguide  is  more 
effective. 

The  maximum  distance  of  sound  propagation  in  the  USC  is  limited  mainly 
by  absorption  in  the  sea  water.  Sound  of  sufficiently  low  frequency,  for  which 


^tp-car 


'to 


(M2) 
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absorption  is  fairly  small,  can  propagate  over  hundreds  and  thousands  of  kilo¬ 
meters. 


1-3.1  Caustics  and  Convergence  Zones 

In  a  ray  diagram,  a  caustic  is  the  envelope  formed  by  the  intersection  of 
adjacent  rays.  It  is  a  region  of  partial  focusing  with  high  intensity.  When  a  caus¬ 
tic  intersects  the  sea  surface,  the  region  at  or  near  the  surface,  within  which  high 
sound  levels  occur,  is  known  as  a  convergence  zone. 

Fig.  1-4  shows  ray  diagrams  for  a  source  at  different  depths  in  an  ocean  with 
the  SSP  shown  in  Fig.  1-3.  These  diagrams  illustrate  the  following  features  for  a 
fairly  typical  ocean: 

•  The  variation  of  relative  insonification  and  shadow  zones  at  different  depths 
and  ranges  for  selected  source  positions. 

•  Xhe  strong  convergence  and  shadow  zones  for  a  source  near  the  surface  or  the 
bottom. 

•  The  importance  of  caustics  in  all  cases. 

As  the  source  approaches  the  USC  axis,  the  width  of  the  shadow  zones 
decreases  and  that  of  the  insonified  zones  increases.  If  the  source  is  on  the  axis 
of  the  USC,  shadow  zones  disappear  at  this  depth.  Several  rays  arrive  at  the 
receiver  if  the  latter  is  not  far  from  the  USC  axis.  Such  multipath  propagation  is 
a  characteristic  feature  of  the  sound  field  in  the  USC. 

The  spatial  and  temporal  (multipath)  structure  of  the  energy  flow  reveals  the 
difficulty  of  locating  an  emitter  in  the  far  field.  As  an  example,  if  the  receiver  is 
located  in  a  shadow  zone,  then  the  source  will  not  be  heard;  therefore,  it  will  not 
be  detected.  Also,  if  the  receiver  is  located  in  an  insonified  zone,  then  it  will 
receive  multiple  signals  from  various  directions  over  a  period  of  time,  causing 
ambiguity  in  determining  the  exact  location  of  the  source.  The  time  spread  of 
the  received  signals  is  a  result  of  the  medium  stratification,  where  different  rays 
(paths)  will  travel  with  different  velocities. 


range  (m) 


FIGURE  1-4.  Ray  diagram  for  source  at  different  depths,  as  shown  on  the  graphs,  for  an 
ocean  with  sound  speed  profile  shown  in  Fig.  1-3.  _ 


CHAPTER  2 


Underwater  Source  Localization 

In  recent  years,  scientists  and  engineers  have  shown  an  increasing  interest  in 
merging  underwater  propagation  models  with  signal  processing  algorithms. 
Here,  we  will  present  three  ways  of  tackling  the  source  localization  problem. 
The  first  two  approaches  are  based  on  meticulous  knowledge  of  the  sound  speed 
structure;  the  third  approach  is  a  novel  one,  and  is  based  on  acoustic  calibration. 
The  general  problem  is  to  infer  the  location  of  the  acoustic  source  from  the  mea¬ 
sured  field  at  a  vertical  receiving  array,  Fig.  2-1. 

2-1  Backpropagation  Method 

The  location  of  an  acoustic  source  can  be  determined  by  backpropagating,  in 
the  same  medium,  the  measured  field  on  the  receiving  array.  This  will  lead  the 
backpropagating  wavefronts  to  converge  at  the  source  location. 


r0  R 


FIGURE  2-1.  Backpropagation  method.  In  ideal  situation  P(r,z)  will  focus  at  r=rt 
and  z=z.. 

The  measured  pressure  field,  P  (r<(  zt;r0,  z) ,  can  be  backpropagated,  toward 
the  unknown  source,  by  convolving  its  conjugate  with  the  spatial  impulse 
response  of  the  medium.  Meticulous  oceanographic  information  is  needed  for 
accurate  representation  of  the  spatial  impulse  response  of  the  medium. 
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The  back  propagated  pressure  field,  from  the  receiving  array  through  a 
medium  with  estimated  spatial  impulse  response,  is 

P(r,z)  =  jp*  (rt,  z,;r0,  z)  G  (r,  z;rQ,  z)  dz  (2-1) 

where  P*  (rf,  zt\r0,  z)  is  the  complex  conjugate  of  the  pressure  field  from  the  tar¬ 
get  at  (rt,zt),  as  measured  at  a  vertical  receiving  array  located  at  range  r0; 

6  (r,  z;r0,  z)  is  the  estimated  spatial  impulse  response  of  the  medium  between 
the  point  (r,  z)  and  a  hydrophone  at  (r0,  z) . 

The  integral  in  Eq.  (2-1)  is  performed  along  the  receiving  array  at  range  r0. 
Ideally,  P(r,z)  will  focus  at  the  exact  location  of  the  source,  r=rt  and  z=zr  The 
quality  of  the  focus  depends  on  the  accuracy  of  the  estimated  Green’s  function,  G  (r,  z) . 

2-2  Matched  Field  Processing  Method 

Bartlett  developed  an  algorithm  for  plane  wave  propagation  in  a  homoge¬ 
neous  medium.  Matched-field  processing  (MFP)  is  a  generalization  of  plane- 
wave  beamforming  [7],  where  the  plane-wave  weighting  vectors,  i.e.,  replicas, 
are  replaced  by  solutions  to  the  wave  equation  for  the  particular  environment  of 
interest. 

Matched  field  processing  has  been  investigated  in  underwater  acoustics  for 
more  than  ten  years.  Matched  field  processors  use  the  field  measurements  as  the 
inputs  and  obtain  the  source  coordinates  as  the  output.  This  processor  is  imple¬ 
mented  by  constructing  forward  solutions  of  the  wave  equation  for  many  test 
source  locations  (the  replica  fields)  and  finding  the  best  “match”  to  the  measured 
data.  This  processor  is  ill-posed  [8];  it  is  highly  unstable  to  such  effects  as 
noise,  errors  in  the  input  parameters  and  the  introduction  of  new  data.  The 
matched-field  problem  itself  is  ill-posed  in  that  it  is  very  sensitive  to  errors  in 
the  model.  This  is  usually  referred  to  as  the  “mismatch”  problem. 

MFP  methods  require  accurate  oceanographic  information  about  the 
medium,  such  as  the  sound  speed  structure  between  the  source  and  the  receiving 
array,  the  attenuation  and  scattering  strength,  the  bathymetry,  and  the  character¬ 
istics  and  structure  of  the  ocean  bottom.  Collecting  and  maintaining  an  updated 
version  of  this  information  over  a  dense  grid  (with  range  increment  about  half 
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the  wavelength)  is  a  very  complicated  task.  Therefore,  inaccurate  oceanographic 
information  will  result  in  false  replica  fields,  and  as  a  consequence  MFP  will 
fail  to  provide  a  unique  solution  to  the  source  location  problem. 

An  ambiguity  function  is  defined  in  terms  of  the  measured  and  replica  fields. 
The  peaks  of  the  ambiguity  function  correspond  to  estimates  of  the  source  loca¬ 
tion  candidates.  Fig.  2-2  shows  a  block  diagram  of  a  conventional  matched  field 
processor.  P  (a)  is  the  measured  pressure  field  from  the  source  along  the  receiv¬ 
ing  array;  P  (a)  are  the  replica  fields  generated  by  a  numerical  wave  propaga¬ 
tion  model  such  as  PE;  each  replica  field  represents  the  pressure  field  along  the 
receiving  array  from  a  hypothetical  test  source  at  range  r'  and  depth  z'.  Differ¬ 
ent  algorithms  use  different  rules  to  find  the  ambiguity  surface  S  (a) .  Bartlett,  or 
the  conventional  estimator,  has  the  ambiguity  surface 

S  (a)  =  P*  ( d)RP  ( a )  (2-2) 

where  the  *  denotes  the  complex  conjugate,  a  =  (r\  z')  is  the  estimated  source 
location,  and  R  is  the  cross  spectral  matrix  of  the  received  field  whose  (i,j)  th 
element  is 


R..  =  E[d.d*]  (2-3) 

where  di  denotes  the  discrete  Fourier  transform  of  the  received  field,  P,(a) ,  on 
hydrophone  i.  Assuming  the  measurement  was  accomplished  over  a  period  of 
time  T,  then  the  discrete  Fourier  transform  is  performed  at  the  frequency  being 
used  for  localization  [9]  and  then  the  expected  value  is  calculated  by  averaging 
over  a  sequence  of  time  frames. 

The  minimum  variance  method  (Capon  estimator)  is  given  by 

S  (a)  =  1/  {ft  (b)R~'P  (3)}  .  (2-4) 

In  case  of  normal  distribution,  the  Capon  estimator  is  usually  referred  to  as  the 
maximum  likelihood  (ML)  estimator[10]. 
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FIGURE  2-2.  Conceptual  diagram  of  the  conventional  matched  field  processor  method. 

Figs.  2-3  through  2-6  represent  the  SSP  and  the  ambiguity  surface  for  a  Bar¬ 
tlett  MFP  for  a  target  located  at  the  axis  of  the  channel  and  at  range  100  km  from 
the  receiving  array.  The  acoustic  frequency  being  used  for  localization  is  100 
Hz.  The  ambiguity  surface  is  plotted  over  50  km  (75  km  to  125  km  from  the 
receiving  array).  The  simulation  was  done  using  a  wide  angle  parabolic  equation 
(PE)  code,  FEPE,  where  the  ocean  was  assumed  to  be  range  dependent  to  simu¬ 
late  the  measured  field.  A  range  independent  SSP  was  used  to  generate  the  rep¬ 
lica  fields.  The  range  dependence  was  modeled  as  a  linear  variation  of  the  sound 
speed  at  the  axis  on  the  order  of  one  meter  per  second  per  hundred  kilometers.  It 
is  evident  from  Fig.  2-4,  where  the  red  is  a  high  intensity  area  and  blue  is  a  low 
intensity,  that  the  focus  is  blurred,  and  the  volume  is  populated  with  false 
images.  Fig.  2-5  shows  a  horizontal  slice  of  the  ambiguity  surface  at  the  source’s 
depth.  High  sidelobes  indicate  the  existence  of  images  at  ranges  other  than  the 
target’s  range,  in  this  case  100  km  from  the  receiving  array.  Fig.  2-6  is  the  same 
as  Fig.  2-5,  except  the  scanning  is  performed  over  a  cylinder  of  radius  5  km 
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FIGURE  2-4.  Color  Image  of  the  ambiguity  surface  for  the  MFP.  The  target  is  at  range 
100  km  and  depth  1000  m.  The  color  bar  shows  the  linear  scale  of  the  normalized  intensity 
grading.  


FIGURE  2-5.  A  horizontal  slice  of  the  MFP  ambiguity  surface  at  z-z,-1000  m. 


2-3  The  Holographic  Approach 

The  mismatch  in  the  MFP  methods  can  be  reduced  by  using  a  reference 
source  near  the  target.  This  approach  brings  the  source  into  focus  by  calibrating 
most  of  the  ocean  between  the  source  and  the  receiving  array.  The  holographic 
method  has  its  genesis  in  optics. 

2-3.1  Holography  in  Optics 

An  outstanding  problem  in  optics  has  been  the  desire  to  minimize  the  effect 
of  the  medium  on  the  resolution  capability  of  imaging  systems,  Goodman,  Hunt- 
ley,  Jackson,  and  Lehman  [1]  has  proposed  the  use  of  a  reference  source  in  the 
vicinity  of  the  object,  so  that  the  object  and  reference  rays  received  at  any  one 
point  on  the  recording  plane  (hologram)  pass  through  identical  portions  of  the 
perturbing  medium.  The  distortions  of  the  reference  and  object  waves  are  thus 
nearly  identical  as  illustrated  in  Fig.  2-7.  The  reference  and  object  waves  reach¬ 
ing  the  hologram  may  thus  be  written  as  PR  (x,  y )  exp  [j\ 3  (x,  y )  ]  and 
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P  t  (x,  y)  exp  [j'P  ( x ,  y)  ] ,  where  exp  UP  (x,  y)  ]  is  the  transmittance  of  the  aberrat¬ 
ing  medium.  The  hologram  records  the  average  intensity  of  the  total  field 

+  +  (2-5) 

where  P*  is  the  complex  conjugate  of  the  object  field,  and  P K  is  the  reference 
field.  The  interference  pattern  between  the  two  is  free  from  the  distorting  effects 
of  the  aberrating  medium.  This  technique  will  work  over  a  limited  object  field, 
since  points  too  far  from  the  reference  may  produce  waves  with  aberrations  dif¬ 
ferent  than  those  of  the  reference  wave. 

A  Hologram  is  a  film  [11]  with  transmittance 

T  =  Kly  (2-5)a 

where  K  is  a  positive  constant,  yis  the  film  contrast  factor1,  and  I  is  the  recorded 
intensity.  Therefore,  coherent  illumination  of  the  hologram  yields 

‘r*  --  +  +  +  (2-6) 

where  we  assumed  y  =  l .  PK  is  the  reference  field  in  the  absence  of  the  aberrat¬ 
ing  medium.  The  fourth  term,  in  Eq.  (2-6),  can  be  neglected  under  the  assump- 
tion  that  the  reference  intensity  is  greater  than  the  object  intensity,  »  |Pt|  . 
The  result  is  a  pair  of  twin  images  of  the  object  symmetrically  located  about  a 
bright  image  of  the  reference. 

Fig.  2-8  illustrates  experimental  results  [1]  obtained  using  conventional  and 
wavefront  reconstruction  imaging  techniques.  The  object  is  a  transparency  con¬ 
taining  both  a  sign  and  a  bright  reference  point,  as  shown  in  Fig.  2-8(a).  The 
aberrating  medium  was  common  shower  glass.  Fig.  2-8(d)  shows  the  wavefront 
reconstruction  images  obtained  in  the  presence  of  the  aberrating  medium.  For 
comparison,  Fig.  2-8(b)  shows  the  conventionally  formed  images,  obtained  by 


1.  7  is  the  slope  of  the  linear  region  of  the  Hurter-Driffield  (H&D)  curve  that  relates  the  photographic  den¬ 
sity  to  the  logarithm  of  exposure.  A  film  with  large  y  is  called  a  high-contrast  film,  while  the  one  with  low  y 
is  called  a  low-contrast  film. 
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replacing  the  film  with  a  positive  lens  which  directly  forms  the  image  in  the 
presence  of  the  aberrating  medium.  Resolution  in  the  conventionally  formed 
image  has  been  totally  destroyed  by  the  aberrations.  The  wavefront-reconstruc¬ 
tion  images,  however,  have  suffered  minor  degradations. 


Random  Medium 


Hologram 

Recording  step===»  1  =  |  (p *  +  p ,)  |2  =  **  +  7, +  p  r  pt  +  prPi 
Reconstruction===»  IP  K  =  I  Rp  r  +  I  Rp ,  + p  rp  Rp* 


FIGURE  2-7.  Conceptual  diagram  of  the  wavefront  reconstruction  (holographic)  imaging 
method  in  optics.  _ _ _ 
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FIGURE  2-8.  The  effect  of  random  media  on  the  resolution  of  imaging  systems;  (a) 
conventional  imaging  with  no  random  medium  present,  (b)  conventional  imaging  in  the 
presence  of  random  medium  between  the  object  and  the  film,  (c)  wavefront  reconstruction 
method  with  no  random  medium  present,  (d)  wavefront  reconstruction  method  with 
random  medium  present  between  the  object  and  the  hologram. 
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2-3.2  Holography  in  the  Ocean 

The  optical  technique  for  wavefront  reconstruction  is  not  directly  applica¬ 
ble  to  long  range  propagation  in  the  ocean,  because  the  sound  speed  profile,  Fig. 
1-3,  varies  significantly  with  depth.  Sound  is  refracted  along  several  paths 
between  the  source  and  the  receiver.  In  the  ocean,  the  received  field  is  the  sum 
of  fields  over  all  paths 


P,  = 

i 


(2-7) 


for  the  reference  source,  and 


p,  =  5/,  o')  <2-8> 

j 

for  the  object,  or  target,  source;  i  and  j  are  the  path  numbers.  The  cross  term  in 
the  holographic  intensity  is 

FV<  =  (2'9) 


FIGURE  2-9.  One  reference  cannot  compensate  for  phase  perturbation  along  all  paths. 


28 


Only  the  cross  terms  for  which  the  object  and  reference  waves  travel  the 
same  paths,  i=j,  will  benefit  from  the  reduction  of  the  phase  perturbation  due  to 
the  random  media.  Since  sound  propagates  to  long  distances  in  the  ocean  over 
many  paths  and  since  each  of  these  paths  is  distorted  in  a  different  way,  a  differ¬ 
ent  phase  correction  must  be  provided  for  each  path.  This  can  be  accomplished 
by  the  use  of  multiple  reference  sources  deployed  at  the  same  range,  but  distrib¬ 
uted  in  depth. 

We  shall  show  that  the  appropriate  holographic  combiner  is  given  by  the 
operation 


H(zr)  =  J  P*  R(z,zR)P,(z)dz  (2-10) 

ho  to  g  ram 

where  z  is  the  position  along  the  receiving  array  (hologram),  and  zRis  the  depth 
of  the  associated  reference  source. 

2-3.3  Wave  Equation  in  a  Range  Dependent  Ocean 

The  acoustic  pressure  P  in  the  ocean  is  a  solution  of  the  Helmholtz  equation, 

2 

9  (0 

VZP  +  —P  =  0  ,  (2-11) 

c 

with  the  appropriate  boundary  conditions,  acoustic  frequency  ox  and  sound 
speed  c.  Consider  the  solution  when  the  depth  is  uniform,  the  medium  is  range 
independent,  the  acoustic  characteristics  of  the  sea  bottom  do  not  change  with 
range,  and  the  sound  speed  is  a  function  of  depth  only,  c=c(z).  The  solution  to 
the  Helmholtz  equation  is  easily  obtained  [12]  by  the  method  of  separation  of 
variables 


P(r,z)  =  tf0(l)  ($r)v  (z)  (2-12) 

where  H 0(1^  (£r)  is  the  Hankel  function  of  the  first  kind  zero  order  (represents 
the  range  dependent  part  of  an  outgoing  wave);  for  large  argument,  %r  »  l ,  it  has 
the  asymptotic  form: 
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«o<0<^>  (2-13) 

where  t,  is  the  separation  constant.  The  depth  function  satisfies  the  wave  equa¬ 
tion 


-4v(z)  +  [*2(z)  -£2]y(2)  =  0 
dz 


and  the  boundary  conditions 

P  (r,  0)  =  0 


(f- P ) 

az 


=  0 


l  =  k 


pressure  release  surface 
rigid  bottom 


(2-14) 


Representing  the  solution  as  the  sum  of  normal  modes,  we  get 

P{r,z )  =  iny£Vm(‘zs)Vn(2)Ho)  ^mr)  (2-15) 

m 

where  \|/m(z)  are  orthonormal  eigenfunctions,  and  are  horizontal  wavenum¬ 
bers  (eigennumbers). 

We  are  concerned  about  the  generalization  of  this  theory  to  the  case  when 
the  characteristics  of  the  ocean  acoustic  waveguide  vary  with  range.  Important 
cases  include: 

•  when  sound  propagates  in  the  coastal  wedge,  where  variation  of  sea  depth  is 
significant, 

•  when  sound  waves  cross  frontal  zones  in  the  ocean,  such  as  currents,  gulf 
stream,  cold  or  hot  eddies,  Kuroshio,  etc., 

•  when  sound  propagates  over  large  ranges  of  the  order  of  thousands  of  kilome¬ 


ters. 
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2-3.3. 1  The  Reference  Waveguide  Method 

To  solve  the  Helmholtz  equation  for  a  range  and  depth  dependent  medium, 
c=c(r,z),  we  will  use  a  reference  waveguide  at  fixed  r=r,  where  c=c(r!,z).  Note 
that  the  depth  depends  on  r,  h  =  h  (r) ,  but  the  reference  waveguide  correspond¬ 
ing  to  the  range  r,  has  a  horizontal  bottom  at  depth  h  =  h(rx) . 

For  the  reference  waveguide  at  an  arbitrary  r,  there  is  a  total  system  of 
orthogonal  eigenfunctions  [13],  \j/m(r,z),  m=l,2,...,  depending  on  z  (r  is  taken  as 
a  parameter)  and  satisfying  the  equation 

2 

(r,  z)  +  [k2  (r,  z)  -  }  y  (r,  z)  =0  (2-16) 

2  **  In  in 

dz 

and  the  boundary  conditions  at  the  surface  and  the  bottom.  Here  ^=^00  are  the 
eigenvalues  of  the  reference  waveguide  at  r  as  illustrated  in  Fig.  2-10.  The  pres¬ 
sure  field  for  a  given  r  can  be  represented  as  an  expansion  in  terms  of  the  eigen¬ 
functions 


P(r,z)  =Xrf(^)Vm(r,z)  (2-17) 

« 

where  Bmare  modal  constants  for  a  range  independent  stratified  waveguide.  Sub¬ 
stituting  P  (r,  z)  into  the  Helmholtz  equation  we  get 


-V  (r,z) 
dz2 


IX".'1’  («.’•> rrv,  <«•.«> 


(2-18) 


V  rP  (r,  ,)  =  0°  (S')  +  Wrlt.H  ^  (S.r)  1  ?,V„  + 

m 


(2-19) 


_  -V  *•  * 

where  V  =  -2—  +  -2—  ,  and  for  brevity  we  omit  the  functional  arguments. 
dx2  dy2 

Taking  advantage  of  the  orthonormal  property  of  the  eigenfunctions  in  each 
reference  waveguide, 
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A 

Jvm(r.z)  *.(»-,  2)  rfz  =  8mH  ,  (2-20) 

o 

substituting  Eq.  (2-18)  and  Eq.  (2-19)  in  Eq.  (2-11),  and  considering  Eq.  (2-16) 
to  substitute  for  and  subsequent  multiplication  by  (r,  z)  followed  by 

integration  over  i^from  0  to  h,  yields  the  following  set  of  coupled  differential 
equations 


(1) 


(2-21) 


2>.«o 


2-3.3.2  Adiabatic  Approximation 

The  right-hand  side  of  Eq.  (2-21)  is  small  if  the  properties  of  the  waveguide 
vary  sufficiently  slowly  with  horizontal  range  r  [14].  If  we  assume  the  right- 
hand  side  is  zero,  then  the  equations  for  the  normal  modes  are  uncoupled.  Each 
mode  propagates  in  the  waveguide  independently  of  the  others.  This  approxima¬ 
tion,  when  the  interaction  between  modes  is  neglected,  is  called  the  adiabatic 
approximation. 

2  If)  r) 

In  cylindrical  coordinates  Vr  =  *  then  the  adiabatic  approxima¬ 

tion  of  the  wave  equation  is 

iA(r±<t>  (r))  +  ^2  (r)0  (r)  =  °  (2-22) 

r  dr  dr  m  m  m 


where  we  use 


<t>m  (r) 


*  A(1)U.  r) 


(2-22)a 


If  we  let 
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Fm(r)  =  rl/2«Dm(r)  ,  <2-22)b 

then  substituting  in  Eq.  (2-22),  and  neglecting  the  terms  with  1/r2  when  solving 
for  the  far  field,  4mr  »  1  >  wc  Set 

=  -en(r)Fm(r)  .  (2-23) 

dr 

If  the  eigennumbers  were  range  independent,  £(r)=£,  then  Eq.  (2-23)  has  the 
solution 


F  (r)  =  Aexp  [±f£r]  (2-24) 

where  A  is  a  constant  of  integration,  to  be  found  from  the  boundary  conditions. 

Dividing  the  range  of  interest  into  intervals  and  assuming  that  £  is  constant 
in  each  interval.  Fig.  2-10,  and  matching  the  solution  at  the  boundaries  of  the 
neighboring  intervals,  then  we  get 


F  (r0,  r) 


^  *5(ro)  ( r~ro ) 


wherer0  <  r  <  r, 


F  (r J,  r) 


where  rx<r<r2  ,  and  by  matching  the  solutions  at  rl5  we  get 

Aj  =  AeU'o)  (ri_ro)  =  A^(,o)A°  , 

and  in  general  for  r{<  r  <  rl  +  x  we  have 


(2-25) 


(2-26) 


= 


A  e 


‘^(ro)  A0  +  *4(ri)A]  +  ...  + 


(2-27) 


33 


ri 


r2 


$2 


FIGURE  2-10.  Dividing  the  range  between  the  source  and  receiver  into  N  incremental 
intervals,  and  the  eigennumber  within  each  interval  is  constant.  Here  we  assumed  that  the 
source  is  at  range  r0.  _ 

Letting  the  range  increments  be  infinitesimally  small,  the  summation  is  equiva¬ 
lent  to  integration,  and  the  range  dependent  function  becomes 


F  (r)  =  AmexP 


(r)  dr 


t-'o 


(2-28) 


Taking  into  account  all  normal  modes,  the  pressure  field  is 


P  (r,  z) 


XA-v-<r’,)r  'n‘xp 


■  r o 


(2-29) 


For  a  stratified,  horizontally  homogeneous  waveguide,  under  the  assumption 
that  %  r  »  l  where  the  asymptotic  form  for  the  Hankel  function  can  be  used,  the 
pressure  field  in  Eq.  (2-15)  becomes 


l/2 

P(r,z)  =  (|)  e-l{n/A)'£vm(z)ym(zi)  ttnr) 


-1/2 


exp  r] 


(2-30) 


At  relatively  small  ranges  r  from  the  source  (but  large  compared  with  the  wave¬ 
length)  the  variation  of  the  medium  properties  in  the  horizontal  does  not  affect 
the  sound  field,  and  ^  can  be  considered  constant.  Thus,  for  these  ranges  Eqs. 
(2-29)  and  (2-30)  must  be  identical.  This  will  yield 
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-  1/2 

2  -i(7t/4)  .,-1/2 

“  (->  C  Vm(0,2,)Sm 


(2-31) 


Then  the  acoustic  pressure  for  a  range  dependent  ocean  is 


P  (r,2) 


(  — ) 
nr 


1/2 


-i  (Jl/4) 


£  (0,  zt)  (r,  z)  (r)  ] 


-1/2 


exp 


-ijtmCr)d~r 
-  ro 


(2-32) 


where  vm  (°> 2,)  is  the  vertical  eigenfunction  in  the  region  close  to  the  source 
and  £m(r)  is  the  eigenvalue  of  the  depth  equation  \|/m(z,r)  for  the  refer¬ 
ence  waveguide  corresponding  to  a  given  r. 

2-3.4  Reconstruction  of  the  scattered  field  at  the  virtual  array 

Without  loss  of  generality  we  adopt  a  cylindrically  symmetric  model  of  the 
ocean.  Fig.  2-11,  where  the  axis  of  symmetry  is  at  r0.  The  sound  speed  varies 
with  depth  and  range,  c(z,r).  The  receiving  array  is  at  coordinate  (Zi,r0)  where  z, 
denotes  the  depth  of  the  individual  hydrophones.  The  unknown  object,  at  posi¬ 
tion  (z„rt),  produces  the  scattered  field  when  illuminated  by  the  reference 
source,  or  when  it  radiates  independently.  The  reference  source  is  shown  at  posi¬ 
tion  (rR,zR).  While  the  following  theory  is  presented  with  the  assumption  that 
rR<rt,  we  note  where  appropriate  the  minor  change  required  for  the  case  rR>rt. 

Synthetic  multiple  reference  sources  can  be  generated  by  periodically  lowering  a 
reference  source  at  the  reference  range  rRand  allowing  it  to  transmit  at  each 
incremental  depth,  Fig.  2-11.  This  way  we  create  a  virtual  array  of  multiple  ref¬ 
erences  distributed  over  the  water  column. 
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r0  surface  rR  rt 


FIGURE  2-11.  A  vertical  slice  of  a  cylindrically  symmetric  ocean  model,  with  a  reference 
source  deployed  near  the  target  to  an  array  of  depths.  _ 


2-3.4.1  Generating  the  Holographic  Image  of  the  Scattered  Field 

A  holographic  image  of  the  scattered  field  from  the  target  can  be  generated 
by  combining  the  measured  pressure  field  from  the  target  with  that  from  each 
and  every  reference  source,  and  then  averaging  over  all  the  hydrophones  of  the 
receiving  array.  We  use  normal  mode  theory  in  a  slowly  varying,  horizontally 
stratified  ocean  to  show  that  this  process  corrects  for  the  variability  of  the 
medium. 

Assuming  that  we  can  separately  identify  the  arrival  of  the  reference  source 
and  target  signals  at  the  receiving  array,  then  the  pressure  field  from  the  target 
on  the  ith  receiving  hydrophone  is 
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and  that  from  a  reference  source  at  depth  zr  is 


2  1/2  -■(/  ) 
p*(ro’2,-)  =  (-)  e  ‘  " 
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exp 


dr 


L  >-0 


(2-34) 


The  real  functions  ym(r,z)  are  the  z-dependent  discrete,  or  ‘propagating’  eigen¬ 
functions  for  the  ‘reference  waveguide  problem’  at  position  r,  and  ^(r)  are  their 
real  eigenvalues.  Although  the  continuous  spectrum,  which  fills  out  the 
orthonormal  basis  set,  can  in  general  make  a  contribution  with  complex  ym, 
these  evanescent  modes  are  significant  only  in  the  immediate  region  near  the 
source,  for  ranges  that  satisfy  the  inequality  kr  <  l ,  where  k  is  the  acoustic 
wavenumber. 

The  accumulated  phase  of  a  mode  propagating  between  two  points  is  given 
by  the  integral  of  the  eigenvalue  over  the  range  between  the  two  points,  as  indi¬ 
cated  in  the  exponents  of  Eqs.  (2-33)  and  (2-34).  The  structure  of  these  equa¬ 
tions  reflects  the  relative  positions  of  the  source  and  receiver.  One  set  of 
eigenfunctions,  \| fm,  in  each  sum  depends  on  the  location  of  the  receivers,  while 
the  other  set  depends  on  the  location  of  the  target  or  reference  source.  * 

The  holographic  image  ‘filter’,  H ,  can  be  constructed  by  combining  the 
pressure  from  the  target  and  the  references 


H  (rt,  Zt;r  zR) 


h 

~  J^i  (ri»  Zi>r0  (rR'ZR'r 

0 


(2-35) 
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Taking  advantage  of  the  orthonormality  of  the  eigenfunctions  of  the  reference 
waveguide  at  r  =  r0,  Eq.  (2-20),  the  double  sum  over  two  different  sets  of 
modes  will  be  reduced  to  a  single  sum,  and  the  modal  representation  of  the  holo¬ 
graphic  image  becomes 


H(2R)  = 


■-Z, - ,  ===—exp 


dr 


(2-36) 


H ,  the  holographic  image,  can  be  viewed  as  a  weighted  version  of  P t  along 
the  virtual  array.  In  contrast  with  the  recording  step  in  optical  holography,  where 
we  record  the  total  intensity  \PT\2,  here  we  record  PR*  ( PT  -  PR) ,  where 
p  =  pR  +  p  t  is  the  total  field  from  the  object  and  the  reference  source. 

H  contains  eigenfunctions  calculated  from  a  reference  waveguide  problem 
at  ranges  rR  and  rt.  Also,  it  has  the  same  phase  as  a  field  propagating  from  a 
source  at  r,and  received  at  rR,  the  reference  locations.  This  field  is 


l/2 
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(2-37) 

If  the  pressure  field  generated  by  the  mth  mode  is  p" ,  and  the  holographic 
image  associated  with  the  mth  mode  is  Hm ,  then 


,(71/4) 
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(2-38) 


Therefore,  we  can  write  P  t,  as  the  weighted  sum  of  the  modes  of  the  holographic 
image  H ,  i.e, 
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Vm(r*>zJ{)  J  H  (rR>  Z*)  Vm  (TJt>  Z*)  dzR 

"  0  (2-39) 


From  this  relation,  we  can  see  that  given  the  holographic  image  at  the  vir¬ 
tual  array,  one  can  calculate  the  pressure  field  from  the  scatterer  as  if  it  was 
measured  at  the  virtual  array. 

The  main  significance  of  H  is  that  it  is  made  up  of  measured  data  where  the 
weighting  of  Pt  by  PR  filters  out  the  effect  of  the  oceanic  variability.  Further¬ 
more,  the  phase  of  each  modal  component  of  H  is  the  accumulated  phase 
between  the  reference  sources  and  the  target.  This  phase  is  independent  of  the 
variability  of  the  index  of  refraction  between  the  source  and  the  receiving  array. 


2-3.4.2  Localization  Algorithm 

Assume  the  ocean  is  frozen  during  the  measurements.  The  localization  algo¬ 
rithm  consists  of  the  following  steps: 

1.  Measure  the  pressure  field  from  the  target  along  the  vertically  distributed  receiving 
array,  P,(z(.). 

2.  Deploy  a  reference  source  at  range  rR  and  depth  zRj  and  allow  it  to  transmit. 

3.  Measure  the  pressure  field  from  the  reference  source.  Its  conjugate  is  P*  (zf,  zR) . 

4.  Multiply  the  pressure  from  the  target  by  the  conjugate  of  the  pressure  from  the  refer¬ 
ence  source  at  each  and  every  hydrophone  of  the  receiving  array. 

5.  Take  the  average  of  the  previous  step  by  adding  the  result  of  multiplication  on  each 
hydrophone  and  dividing  by  the  number  of  hydrophones  in  the  receiving  array  giving 

6.  Deploy  the  reference  source  to  a  new  depth  and  repeat  step  2  through  5. 

7.  When  all  reference  depths  are  used  then  we  have  the  total  holographic  image, 

H  ( zR  j ,  zR  2, ... ,  zR  N ) ,  of  the  target  as  if  it  was  recorded  on  the  virtual  array. 

8.  While  deploying  the  reference  source,  we  measure  the  SSP  of  the  ocean  at  the  virtual 
array. 


9.  Assume  the  ocean  is  range  independent  between  the  virtual  array  and  the  target,  and  it 
has  the  SSP  obtained  in  step  8. 

10. Using  a  propagation  model,  such  as  wide  angle  PE,  compute  the  pressure  field  on  the 
virtual  array  from  a  hypothetical  point  source  at  position  (r',z')  in  the  vicinity  of  the 
reference  and  target,  giving  G  ( rR ,  zR  ;r',  z') ,  the  spatial  impulse  response  between 
the  hypothetical  location  and  the  virtual  array. 

11.  Multiply  the  pressure  from  the  hypothetical  position,  G,  by  the  conjugate  of  the  holo¬ 
gram  H,  and  average  over  all  elements  of  the  virtual  array  giving  the  reconstructed 
field. 


hVA 

Q(r',z')  =  j  G  (rg,  zR;r\  z')  H*  (rK,  zR-,rt,  zt)  dzK  (2-40) 

o 

12. Repeat  the  steps  10  through  11  exhausting  all  possible  candidate  locations  of  the  target, 
yielding  a  matrix  of  Q  in  range  and  depth. 

13.  The  source  location  is  where  the  magnitude  of  Q  attains  a  maximum,  then 

(r\  z ')  =>  (r t,  zt )  . 


The  above  steps  can  be  summarized  in  a  block  diagram  format,  as  shown  in 
Figs.  2-12  and  2-13,  where  Fig.  2-12  shows  a  simplified  block  diagram  for  the 
holographic  array  processing  algorithm,  and  Fig.  2-13  shows  a  more  detailed 
block  diagram  for  the  HAP  algorithm.  In  Fig.  2-13  we  assumed  N  receivers  and 
M  reference  sources.  The  STD  is  an  instrument  that  measures  salinity,  tempera¬ 
ture,  and  depth;  this  data  can  be  used  to  calculate  the  sound  speed  profile  at  the 
virtual  array,  as  described  in  chapter  1. 

The  HAP  method  assumes  no  communication  link  between  the  reference  and 
the  receiver,  and  the  data  processing  takes  place  at  the  receiving  site.  In  this 
case,  the  collected  STD  data  will  be  useful  only  for  updating  the  oceanographic 
data-base  for  that  region  for  future  implementation. 


o  o  o** •  o 


Recording  Step 


FIGURE  2-13.  Detailed  block  diagram  for  the  HAP  algorithm. 


•  •  • 
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2-3.4.3  Focusing  Criterion 

In  our  simulation  we  assumed  that  the  ocean  is  range  independent  between 
the  target  and  the  reference  array.  The  estimated  eigennumbers  and  eigen¬ 
functions  \\fm  are  those  of  the  reference  waveguide  at  the  virtual  array.  Then  the 
pressure  field  from  a  test  source  at  location  ( r\z '),  as  calculated  at  the  virtual 
array,  is 


,  IN  ,2  i 


dr 


.  (2-41) 


Here  we  assumed  that  rR  <  r',  although  a  change  in  sign  in  the  exponential  would 
describe  the  case  where  rR>  r'. 

Since  H  is  an  approximation  to  the  field  from  the  target  on  the  virtual  array, 
and  G  constitutes  an  approximation  to  the  spatial  impulse  response  (Green’s 
function)  of  the  medium  between  the  target  and  the  virtual  array,  backpropagat- 
ing  H  through  this  medium,  it  focuses  on  the  location  of  the  target.  Representing 
Q  in  the  normal  mode  formalism  we  get, 


Q(r',z’)  = 


1  'f,n(rr 

'  r,  r 

exp  I  /J\(r)  dr-ijim(r)  dr 


(2-42) 


A  normalized  version  of  Q  is 


Qn (r'-z')  =  JrRr:\r'~rR\  (^)  2(r’’z’) 


-3/2 


(2-43) 


By  assuming  that  the  ocean  is  range  independent  on  the  scale  of  the  separation 
between  the  reference  and  the  target,  in  the  region  of  the  source,  we  can  write 
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c  (2,  rR)  -  c  (2,  r')  -  c  (2,  r().  Hence  (r',  2)  «=  \f  m  (rt,  z)  »  (rR>  z) ,  permitting 
the  weak  approximation 


lim  v„  (»>*,)  2')  -  Jim  V.  (»•*.  2,)  (»V  *') 

X  — »x,  x  — >  X, 


.  (2-44) 


In  the  limit  2'  -»  2,,  QN(r',z')  reaches  a  maximal  value,  since  it  starts  as  a 
sum  of  terms  of  varying  sign  and  becomes  a  sum  of  positive  numbers.  By  identi¬ 
fying  z'  and  r'  which  maximize  QN,  we  have  located  a  candidate  position  of  the 
unknown  source,  without  knowing  the  oceanic  variability  between  the  receiving 
and  virtual  arrays. 

This  algorithm  eliminates  the  need  to  simulate  the  index  of  refraction, 
except  in  the  region  of  the  reference  array.  Ocean  simulation  is  reduced  from 
basin  scales  to  scales  characteristic  of  the  neighborhood  of  the  sources. 


CHAPTER  3 


Performance  Analysis  of  the  Holographic  Array  Processor 

In  the  following  we  provide  a  performance  analysis  of  the  Holographic 
Array  Processing  (HAP)  algorithm  in  a  range  dependent  ocean.  For  this  purpose, 
we  consider  three  models  of  the  ocean: 

•  A  hypothetical  ocean  with  range  dependent  bP-bilinear  Sound  Speed  Profiles  (SSPs). 

•  A  hypothetical  ocean  with  range  dependent  exponentially  stratified  (Munk)  SSPs. 

•  A  realistic  ocean  with  SSPs  from  the  North  Pacific  experimental  data. 


3-1  Theoretical  Analysis 

This  analysis  is  done  for  a  range  dependent  ocean  with  N2-bilinear  SSP 
which  is  analytically  attractive  and  representative  of  a  wide  variety  of  regions  in 
the  ocean.  If  the  estimated  horizontal  wavenumber  |m(r)  is  very  close  to  the 
actual  Z,m(r)  in  the  scanning  region  (between  the  target  and  the  virtual  array), 
then  the  phase  term  in  the  expression  of  QN 

r,  f 

t1„  =  0-1) 

rR  rR 

will  be  close  to  zero  for  all  modes  as  r'  approaches  rt,  and  the  summation  in  Eq. 
(2-42)  includes  only  positive  quantities.  Then,  QN  will  attain  a  local  maximum, 
identifying  a  candidate  for  the  location  of  the  unknown  source.  Uncertainty  of 
the  estimated  parameters,  for  the  ocean  between  the  reference  array  and  the  tar¬ 
get,  will  lead  to  source  localization  errors.  The  phase,  which  is  a  function  of  the 
horizontal  wavenumber,  provides  information  about  the  focus  in  range,  while 
the  amplitude  (vertical  mode  functions)  provides  information  about  the  focus  in 
depth.  The  following  analysis  will  be  centered  on  the  focus  error  in  range. 
Therefore,  we  need  explicit  expressions  for  the  estimated  and  exact  eigennum- 
bers,  £m(r)  and  \m(r)  *  pertaining  to  the  scanning  region  between  the  target  and 
reference  source. 
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For  a  range-dependent  ocean,  the  index  of  refraction  varies  with  distance 
and  depth.  We  normalize  the  index  of  refraction  by  a  reference  sound  speed.  For 
convenience,  the  reference  sound  speed  is  assumed  to  be  ca0,  the  sound  speed  at 
range  r=0  and  depth  z=za,  where  za  is  the  depth  of  the  channel  axis.  Thus 

w<^)  =  cfe)  <3-2) 

is  the  normalized  index  of  refraction.  Also,  we  assume  that  the  sound  speed  at 
the  axis  varies  linearly  with  range, 

ca(r)  =  Ca0  +  Sr  (3-3) 

where  S  is  the  rate  of  variation,  and  is  measured  in  [(m/sec)/100  km]. 

Consider  the  following  model  for  the  index  of  refraction  N  (see  Fig.  3-1): 

•  The  vertex  of  the  profile  is  fixed  at  depth  za  for  all  ranges, 

•  the  depth  of  the  ocean  zx  is  constant, 

•  the  sound  speed  at  the  surface  and  the  bottom  of  the  ocean  is  constant  and  is  cu  and  q 
respectively, 

•  we  assume  that  cu=  q, 

•  the  N2-bilinear  SSP  is 


N2(z,r)  =  N2(za,r)  +Su{l0)(z-za) 


(3-4) 


where  is  the  slope  of  the  upper  or  lower  portion  of  the  N2-bilinear  profile, 


S„(r)  =  -- 
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(3-5) 


We  use  the  WKB  approximation,  which  is  valid  in  the  geometric-acoustic 
limit  [15],  to  calculate  an  approximation  to  the  eigennumbers  %  .  Following  the 
procedure  presented  in  appendix  I,  the  vertical  wavenumbers  have  a  phase  accu¬ 
mulation  of 


8(z.zwp)  = 


—  f 


N(z,r)  2  a 

(ttt - r)  -  COS>  dz 

yN(z,r)' 


utp 


(3-7) 


resulting  from  the  upward  and  downward  traveling  waves.  Since  these  waves  are 
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out  of  phase  by  nil  radians  upon  reflection  from  the  turning  points,  the  reso¬ 
nance  condition  that  isolates  the  eigenmodes  is 

6  +  <3*8) 

where  zmp  and  zUp  are  the  upper  and  lower  turning  points  respectively  for  each 
mode.  Their  dependence  on  the  mode  number  (m)  is  suppressed  in  this  notation. 
These  turning  points  can  be  found  by  setting  the  integrand  in  Eq.  (3-7)  to  zero 
which  is  equivalent  to  defining  the  turning  points  where  the  vertical  wavenum- 
ber  (Ym)  vanishes;  see  Fig.  3-2.  For  the  N2-bilinear  profile,  we  substitute  Eq.  (3- 
4)  in  the  integrand  of  Eq.  (3-7)  and  equate  to  zero.  After  some  mathematical 
manipulation,  we  get 

zut  =  za~  zaN2  (za’  r )  C  W  si°2  1  (3*9) 

zlt  =  2 a  +  <z/  ~  n1  (za’  r)  C  (r)  sin2  1 1 (r)  1  (3*10) 

where 


Y=0 

FIGURE  3-2.  Variation  of  the  vertical  wavenumber  between  the  turning  points. 
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The  range  dependent  mode  propagation  wavenumber,  £m,  can  be  written  as 

Mr>  =  ^C°S(p-(r)  (3*12) 


where  <pm(r)  is  the  range  dependent  ray  exit  angle  associated  with  the  mode  m. 
Substituting  Eq.  (3-8)  into  Eq.  (3-7),  we  get 


(m  +  -)7c 


z « 


(3-13) 


depth  (m) 

FIGURE  3-3.  The  upper  (solid  line)  and  lower  (dotted  line)  turning  point  depths  as  a 
function  of  mode  number. 
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After  some  mathematical  manipulation,  we  find  coscp^  (r)  as  follows: 


cosz(pm  (r)  =  1- 


cj') 


-,2/3 


z/CoNz(zfl,r)C(r)J 


(3-14) 


and  let 


a  (m) 


rl(m+\)n 

ZyCO 


-.2/3 


(3-15) 


P(r) 


cfl(0 


-,2/3 


7V2(zfl,  r)C(r) 


(3-16) 


Therefore,  we  can  write 


5m(r)  =  FIT)  [1-“(m,P(r)]1/2-  (3-17) 

An  excellent  approximation  to  Eq.  (3-17)  can  be  derived  by  noting  that 
Snell’s  law  determines  the  upper  limit  on  the  ray  exit  angle  <pm  and,  therefore, 
on  the  number  of  propagating  modes 


Q 

cos  [<p"“  (r)  ]  =  =  N  (z„  r)  -  1.0.  (3-18) 

Since  N  is  close  to  one,  <p™ax  is  small,  on  the  order  of  13°.  Applying  Eq.  (3-18) 
to  Eq.  (3-14)  and  Eq.  (3-17)  we  conclude  that 


a  ( m )  (3  (r)  «  1  . 


(3-19) 
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We  can  expand  Eq.  (3-17)  using  the  binomial  series,  producing 

^(r)  =  • 


(3-20) 


This  approximation  is  nearly  exact  (within  106)  for  the  ocean  parameters  we 
have  considered.  A  linearized  version  of  £m(r)  in  the  vicinity  of  the  virtual 
array  is 


£  (r) 


(3-21) 


0) 


Evaluating  these  terms,  and  noting  that  ka(r)  —  — - r-,  we  obtain  an  analyti- 

cal  expression  for  £  (r) . 


Ic  ( ^ 


1  +a(m) 


1B,  ,  2fP<'>> 
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(r~rR) 


(3-22) 


Fig.  3-4  compares  the  WKB  expression  for  the  eigennumber,  Eq.  (3-17), 
with  its  approximation,  Eq.  (3-22),  for  different  modes.  From  the  graphs  it  is 
noticeable  that  the  difference  between  the  two  expression  is  at  its  maximum  for 
the  highest  mode  number  at  the  largest  distance,  and  this  is  of  the  order  lO^m1. 

In  practice,  one  can  measure  N(z,  rR )  using  a  conductivity-temperature- 
depth  (CTD)  casts  near  the  reference  sources.  This  information  is  used  in  a  wave 
propagation  model  to  calculate  QN.  For  our  analysis,  these  data  are  used  to  esti¬ 
mate  the  mode  propagation  wavenumber,  |m(r) ,  between  the  target  and  the  vir¬ 
tual  array, 


I mOO-$m(»>)  • 


(3-23) 


When  numerically  evaluating  the  focusing  condition,  Eq.  (3-1),  we  use  the 
WKB  solution,  Eq.  (3-17),  for  %  .  The  same  expression  is  evaluated  at  rRfor  the 
estimate  |m.  When  analytically  evaluating  Eq.  (3-1)  we  use  the  linearized 
expression,  Eq.  (3-22),  for  t,m  and  the  same  expression  evaluated  at  r  =  rR  for 
the  estimated  Evaluating  the  integrals  in  Eq.  (3-1),  we  define 
A  =  r'  -rR,  =  rt~rR,  8£  =  A-A^  Then  we  can  transform  the  condition 
■q  =  0  into  an  expression  for  the  focusing  error  (in  range)  of  the  unknown 
source  as  a  function  of  mode  number, 


6£(m)  = 


ka  > 


(1  +a(m) 


c(za,rR) 


(3-24) 


3.1.1  Discussion  of  the  Localization  Error  for  a  Range  of  Physical  Conditions 


Figs.  3-5  through  3-8  show  the  dependence  of  the  focusing  error  5e  on  the 
mode  number  for  an  ocean  with  the  following  parameters  (see  Fig.  3-1):  Zj=5 
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km,  c,o=1500  m/sec  and  cu=clo=1540  m/sec,  frequency  100  Hz,  a  variation  of  tar¬ 
get/reference  source  separations,  and  a  range  of  horizontal  gradient  in  the  sound 
speed  at  the  sound  channel  axis.  In  all  of  these  figures,  we  considered  contribu¬ 
tions  from  propagating  modes  only. 

In  all  figures  where  values  of  the  environmental  parameters  are  chosen  to 
represent  typical  range  variations  encountered  in  the  ocean,  the  error  SE(r)  in  the 
estimated  range  of  the  unknown  source  is  nearly  constant  with  mode  number.  In 
particular,  the  variation  in  8E  as  a  function  of  mode  number  is  a  small  fraction  of 
the  wavelength  of  the  acoustic  wave.  The  entire  phase  term  T]m  would  only  vary 
a  fraction  of  a  radian  in  these  cases,  and  therefore,  there  would  be  a  local  maxi¬ 
mum  in  the  localization  estimator  Q  for  the  predicted  value  of  range. 

It  is  also  true  that,  with  oceanographically  motivated  environmental  parame¬ 
ters,  the  error  in  the  predicted  range  of  the  unknown  source  is  a  fraction  of  a 
wavelength.  For  example,  Fig.  3-5  shows  8E  for  At  (the  distance  between  the  ref¬ 
erence  array  and  the  unknown  source)  on  the  order  of  10  and  20  km.  The  value 
of  the  horizontal  gradient  in  sound  speed  is  lO^sec*1.  The  error  quadruples  as  the 
distance  between  the  source  and  the  reference  is  doubled.  However,  the  actual 
values  of  the  error  in  range  is  still  less  than  ten  percent  of  the  15  m  wavelength 
of  the  100  Hz  sound,  and  the  same  is  true  for  the  variation  of  8e  with  mode  num¬ 
ber. 

Figs.  3-6  and  3-7  show  8E  for  target-reference  separation  (A,)  of  10  and  20 
km  respectively,  and  a  range  of  S,  the  horizontal  gradient  in  sound  speed,  from 
the  oceanographically  motivated  value  of  10'5to  four  times  that  value.  Both 
cases  show  that  the  absolute  range  error  and  the  variation  in  error  with  mode 
number  are  small  compared  with  the  acoustic  wavelength.  However,  with  an 
increase  in  S  beyond  4  x  10  the  variation  in  range-localization  error  with  the 
mode  number  becomes  an  appreciable  fraction  of  a  wavelength.  With  this  much 
variation,  the  phase  of  Q  may  not  be  small  enough  for  a  substantial  number  of 
the  propagating  modes  to  add  coherently.  Therefore,  the  argument  for  a  local 
maximum  in  Q  may  not  apply  and  the  algorithm  may  fail.  Further  investigation 
needs  to  be  done  to  determine  the  effect  of  rapidly  varying  horizontal  sound 
speed  gradients. 
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Fig.  3-8  shows  both  numerical  and  analytical,  Eq.  (3-24),  expressions  for  8e; 
the  target-reference  separation  is  20  km,  and  the  series  of  the  horizontal  sound 
speed  gradients  is  the  same  as  in  the  previous  plots.  The  difference  between  the 
analytical  and  numerical  solutions  is  negligible. 

In  Fig.  3-9  we  show  the  frequency  dependence  of  the  range-localization 
error,  8E,  for  target-reference  separation  of  20  km  and  horizontal  sound  speed 
gradient  of  one  meter  per  second  per  100  km.  The  error  approaches  zero  as  fre¬ 
quency  decreases  for  a  fixed  mode  number,  because  the  modes  cover  a  large 
fraction  of  a  vertical  slice  of  the  water  column  as  frequency  decreases.  In  the 
waveguide  that  we  are  considering,  the  largest  change  in  sound  speed  with  range 
occurs  near  the  axis  of  the  sound  channel,  thereby  affecting  the  higher-fre¬ 
quency  localization  procedure  more  strongly  for  a  fixed  mode  number. 

Fig.  3-10  shows  the  phase  term  Tim  for  f=100  Hz,  \=20  km,  and  S=10  5  sec  1 
(a  variation  in  the  axial  sound  speed  of  1  meter  per  second  per  100  km).  The 
numerical  version  of  the  integrand  is  used  for  these  calculations.  The  values  of 
8e  shown  on  the  graph  translate  into  different  constant  values  for  r'  to  be  put 
into  Eq.  (3-  1)  for  rim.  Note  that  Tjm  varies  only  a  fraction  of  a  radian  from  0  for 
the  different  modes.  Therefore,  there  would  be  a  local  maximum  in  the  localiza¬ 
tion  estimator  Qjq,  which  would  indicate  the  predicted  range. 
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FIGURE  3-5.  Focusing  error  8e  (i.e.,  the  difference  between  the  predicted  and  the  actual 
range  of  the  unknown  source)  as  a  function  of  the  mode  number  for  two  values  of  At,  the 
distance  between  the  unknown  source  and  the  reference  array.  _ _ _ 


FIGURE  3-6.  As  in  Fig.  3-5,  but  for  At=10  km  and  a  series  of  values  of  s,  the  horizontal 
gradients  in  the  channel-axis  sound  speed.  _ _ _ _ _ 
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FIGURE  3-7.  As  in  Fig.  3-6,  but  for  A,=20  km. 


FIGURE  3-8.  A  comparison  of  the  numerical  integration  (solid  line)  of  the  localization 
condition  Tim  =  0  with  the  analytical  solution,  Eq.  (3-24).  _ 
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3-2  Numerical  Simulation 

Numerical  simulations  have  been  executed  with  range-dependent,  analytical 
and  measured,  sound  speed  profiles.  An  adiabatic  normal  mode  propagation 
model  has  been  developed,  and  a  wide  angle  wave  propagation  model,  FEPE, 
has  been  modified  for  testing  the  holographic  array  processing  algorithm.  Both 
models  produce  substantial  focusing  at  the  location  of  the  target  source.  The  two 
methods  reveal  slight  differences  for  the  test  positions  (r\z')  away  from  the  tar¬ 
get  location  (r^zj.  This  is  mainly  a  result  of  the  different  approximation  meth¬ 
ods  used  to  solve  the  wave  equation.  While  the  adiabatic  mode  propagation 
model  is  numerically  efficient  for  a  slowly  varying  ocean,  the  PE  model  can 
handle  more  severe  changes  in  the  medium  and  account  for  mode  coupling  terms 
that  are  ignored  by  the  adiabatic  mode  model. 

3-2.1  Exponentially  Stratified  Sound  Speed  Profile 

Munk  [16]  [17]  developed  a  widely  used  mathematical  expression  for  the 
SSP  in  the  ocean. 


c(z)  =  c(za)  [1  +  e(v  +  e  v—  1)]  (3-25) 

where  v  =  —  (z-z  ) ,  z  is  the  depth  of  the  channel  axis,  Hc  is  the  characteris- 
tic  scale  depth,  and  e  is  a  dimensionless  shape  parameter.  Three  parameters  (z„ 
e,  and  v)  control  the  behavior  of  the  Munk  profile.  Fig.  3-11,  shows  a  typical 
exponentially  stratified,  Munk,  SSP,  where  the  axis  depth  is  taken  at  1000  m,  Hc 
is  1000  and  e  is  0.057. 

The  range  dependency  is  introduced  by  assuming  that  the  sound  speed  at  the 
axis  varies  linearly  with  range, 

c  (z,  r)  =  c  (za,  r0)  +  S  ■  r  .  (3-26) 

The  profiles  shown.  Fig.  3-11,  are  separated  by  500  km  with  S=10'5  sec1  which 
corresponds  to  a  change  in  the  axial  sound  speed  of  1  m/sec  per  100  km.  The 
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The  depth-step  at  the  receiving  and  virtual  arrays  is  such  that  there  is  ade¬ 
quate  sampling  for  all  modes,  especially  the  high  order  modes.  The  highest 
order  mode  has  the  largest  spatial  frequency  and,  therefore,  the  smallest  spatial 
wavelength.  Thus,  sufficient  sampling  of  this  mode  will  guarantee  adequate 
sampling  for  the  other  (lower  order)  modes.  The  largest  vertical  wavenumber  is 

fiax  =  ^  sin  6max  (3-27) 

a 

where  is  the  acoustic  wavelength  and  Qmax  is  the  maximum  exit  angle  (corre¬ 
sponds  to  the  highest  order  mode).  The  vertical  wavenumber  can  be  written  in 
terms  of  the  spatial  wave  length  (A,  )  as 

y nax 


Substituting  Eq.  (3-27)  in  Eq.  (3-28)  we  get 
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X 


min 

s 


X 


sin  0' 


max 


(3-29) 


Therefore,  the  highest  order  mode  will  be  adequately  sampled  if  the  depth-step 
is  less  than  or  equal  to  half  the  minimum  spatial  wavelength,  Eq.  (3-29). 

Figs.  3-13  and  3-14  show  the  normalized  magnitude  square  of  Q  as  a  func¬ 
tion  of  the  scanning-source  range.  In  the  case  of  the  PE  model,  Figs,  (b),  we 
assumed  an  absorbing  bottom  and  an  acoustic  frequency  of  50  Hz,  and  a  mesh 
setup  with  increments  of  100  m  in  range  and  50  m  in  depth.  The  same  parame¬ 
ters  for  the  adiabatic  model  were  used,  except  no  bottom  model  is  included, 
since,  only  the  water  borne  modes  are  considered.  In  both  cases  we  assumed  that 
the  target/reference  separation  is  25  km,  so  the  test  source  will  scan  an  area  of 
50  km  with  the  target  between  the  reference  array  and  the  farthest  searching 
point.  On  the  plots  the  virtual  array  is  located  at  range  0  km,  the  target  at  25  km, 
and  the  receiving  array  (not  shown)  at  -450  km. 

A  sharp  focus  of  acoustic  signal  is  observed  at  the  correct  location  of  25  km 
in  Figs.  3-13  and  3-14.  The  width  of  the  central  lobe  is  about  500  m  at  the  1/2 
power  level.  The  highest  sidelobes  are  more  than  10  dB  down  from  the  peak. 
Simulation  for  the  focusing  in  depth  was  done;  see  Fig.  3-15.  This  plot  shows  a 
color  image  of  the  normalized  Q  magnitude,  where  the  color  scale  is  linear  with 
red  representing  high  intensity  and  blue  representing  low  intensity.  This  image 
was  obtained  for  the  same  parameters  as  for  Figs.  3-13  and  3-14,  using  the  PE 
propagation  model  (FEPE).  The  focus  is  clearly  observed  at  the  exact  range  and 
depth  of  the  target.  A  very  sharp  focus  in  depth  at  the  location  of  the  target  was 
found.  The  explanation  for  this  behavior  is  as  follows.  When  we  are  at  the  target 
range,  the  phase  of  each  mode  is  zero,  and  what  remains  in  the  expression  for  Q 
are  the  eigenfunctions,  which  are  oscillatory  functions.  When  we  are  at  the  right 
depth  and  the  right  range,  (r',  z')  =>  (r,,  zt) ,  as  it  was  explained  in  the  previous 
chapter,  the  terms  in  the  summation  are  positive.  When  z'  is  different  from  zt, 
the  summation  tends  to  zero. 
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FIGURE  3-15.  Color  image  of  the  reconstructed  field  using  HAP  with  a  range  dependent  Munk  SSP. 
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3-2.2  Measured  Sound  Speed  Profile 

Real,  measured  SSP  can  provide  a  more  challenging  test  of  the  HAP  method 
than  the  Munk  profile.  Fig.  3- 16(a)  shows  the  region  of  the  North  Pacific  where 
the  experiment  SLICE  89  was  conducted.  Real  data  where  collected  along  the 
three  legs  of  the  triangle  shown  on  the  picture  [  1 8]  [  19].  One  set  of  data  was  col¬ 
lect  during  deployment  in  May  89,  and  another  set  of  data  was  collected  during 
recovery  in  September  89.  This  will  give  six  sets  of  data  with  geographical  and 
seasonal  variations.  In  testing  the  holographic  processor  we  used  the  data  col¬ 
lected  in  May,  along  leg  1-2  (North-South).  The  length  of  this  leg  is  950  km  and 
a  profile  was  taken  every  10  km. 

Three  representative  profiles  are  shown  in  Fig.  3- 16(b),  taken  from  the  May 
data  set:  The  beginning  of  the  path  “1”  (the  solid  line),  the  middle  of  the  path 
(the  dashed  line),  and  500  km  south  of  point  “1”  (the  doted  line).  We  have  a  total 
of  50  profiles.  We  used  FEPE  with  these  profiles  and  assumed  a  target  located  at 
475  km  south  of  point  “1”,  and  a  virtual  array  at  450  km  south  of  point  “1”,  as 
previously  for  the  Munk  profile.  Fig.  3-17  shows  strong  focusing  of  the  HAP, 
similar  to  what  we  have  with  the  Munk  profile.  The  1/2  power  width  is  still  less 
than  1  km,  and  the  sidelobes  are  more  than  10  dB  below  the  focal  point  level. 
The  focal  region  is  strongly  asymmetric,  reflecting  the  strong  irregular  range 
dependency  of  the  SSP. 
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FIGURE  3-16.  (a)  SLICE  89  experiment  course  (triangle),  (b)  Representative  SSP  along 
the  North  -  South  leg  from  the  May  data,  over  a  range  of  500  km:  the  solid  line  represents 
the  SSP  at  point  1,  the  dashed  line  represents  the  SSP  at  250  km  south  of  point  1,  and  the 
doted  line  represents  the  SSP  at  500  km  south  of  point  1. 


CHAPTER  4 


TVuncated  Array 

In  deep  water,  the  ideal  situation  of  having  a  receiving  array  and  reference 
source  that  span  the  water  column  may  not  be  practical.  We  may  have  to  deal 
with  a  receiving  array  that  extends  from  the  surface,  or  just  below  the  surface,  to 
a  certain  depth,  usually  much  shallower  than  the  depth  of  the  ocean.  Similarly, 
we  may  have  to  deal  with  a  reference  source  that  cannot  be  deployed  very  deep. 

In  this  chapter  we  analyze  how  truncating  both  the  receiving  and  the  refer¬ 
ence  source  array  affects  the  performance  of  the  holographic  array  processing 
algorithm. 

4-1  Theoretical  Analyses 

Assume  that  the  receiving  array  extends  from  the  surface  to  a  certain  depth, 
h,  where  h  <h,  and  h  is  the  total  depth  of  the  water  column,  as  illustrated  in  Fig. 
4-1.  Also,  we  assume  that  the  reference  source  spans  the  water  column  from  the 
surface  to  a  certain  depth,  Hr,  where  hR<hR  and  hR  is  the  total  depth  of  the 
water  column  where  the  reference  source  is  deployed. 


FIGURE  4-1.  Sketch  of  the  ocean  model  with  truncated  receiving  and  virtual  arrays. 
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Using  these  definitions,  the  holographic  image  at  the  virtual  array  is 

h 

H  (zR)  =  jPt  ( z\rv  zt)  PR*  ( z;rR ,  zR)  dz  (4-1) 

0 

where  P t  ( z\r{ ,  zt)  is  the  pressure  at  the  receiving  array  from  the  target  located 
at  distance  rt  and  depth  zt,  and  PR  ( z;rR ,  zR)  is  the  complex  conjugate  of  the 
received  pressure  from  a  reference  source  at  range  rR  and  depth  zR.  The  holo¬ 
graphic  image  can  be  rewritten  in  terms  of  the  ocean  depth  h  as, 

h  h 

H  (zR)  =  \Pt  (z;rr  zt)  PR  (z;rR,  zR)  dz  -  (z;r{,  z,)  PR  ( z;rR, ,  zR)  dz  .  (4-2) 

0  h 

The  first  term  represents  the  holographic  image  H  ( rR ,  zR) ,  in  the  case  of  a  full 
length  receiving  array  that  spans  the  water  column  from  surface  to  bottom, 
h  -  h.  Thus,  Eq.  (4-2)  can  be  written  as 

h 

H  {zR)  =  H  (r R,  zR)  -J>,  ( z;rt ,  z()  PR  ( z;rR ,  zR)  dz  .  (4-3) 

h 

The  reconstructed  image,  as  a  result  of  weighting  the  holographic  image  by 
the  pressure  from  a  test  point  source  at  range  r'  and  depth  z',  P  ( rR ,  zR,r\  z'), 
and  averaging  over  the  truncated  virtual  array,  is 


hjt  * 

Q(r',z')  =  fP(rR,zR;r',z')H  (r R,  zR)  dzR  (4-4) 

0 

* 

where  hR  is  the  truncated  depth  of  the  virtual  array.  In  terms  of  the  ocean  depth 
at  the  reference  location,  hR,  the  reconstructed  image  is 


69 


P  (rR,  zR',r',  z')H*  (zR)dzR-  j  P(rR,zR;r',z')H*  ( zR)dzR  .  (4-5) 

ftR 

Substituting  Eq.  (4-3)  into  Eq.  (4-5)  we  get 

^R 

Q  (r\  2 ')  =  J  P  (rR>  ZR’r  >Z  (rR>  ZR )  dzR 

0 

^ r  h 

—  J  dzR^dzP  ( r R,  zR\r  >  z  )  P  t  (2!^»  zt)  Pr  (z’rR>  zr ) 

£  h 

hR 

~\P(.rR,zR,r\  Z')H*  (zR)dzR 

hR 

K  h 

+  J  dzRjdzP  (r R,  ZR-J  > 2  )  P i  (2>f"j>  zi )  P r  (2>*"ft>  zr ) 
hR  h  (4-6) 

The  first  term  in  Eq.  (4-6)  is  the  original  definition  of  the  reconstructed  image 
when  neither  the  receiving  array  nor  the  virtual  array  is  truncated,  and  it  is  equal 
to  Q  (r\  z ') .  The  second  term  gives  the  effect  of  truncating  the  receiving  array, 
the  third  term  gives  the  effect  of  truncating  the  virtual  array,  and  the  fourth  term 
gives  the  effect  of  truncating  both  the  receiving  and  virtual  arrays. 

Letting  h  =  h,  thus  truncating  the  virtual  array  only,  the  second  and  fourth 
terms  are  zero,  and  the  reconstructed  image  will  be 

^R 

Q  (r',  z)  =  Q  (r',  z')  -  jP  (r R,  zR\r\  z')  ff  ( zR )  dzR  .  (4-7) 

tlR 

While,  letting  Hr  =  hR  ,  thus  truncating  the  receiving  array  only,  the  third  and 
fourth  terms  will  be  zero,  and  the  reconstructed  image  will  be 


Q(r',zf)  =  J 
0 
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A r  A 

Q  ( r z')  =  Q  (r',  z')  -  J  dzRjdzP  ( rR ,  zR,r\  z')  P*  ( z;r{ ,  z,)  PR  ( z;rR z*)  .  (4-8) 

0  A 

The  second  term  in  Eq.  (4-7)  and  Eq.  (4-8)  represents  the  disturbance  (error) 
factor  due  to  the  two  truncation  scenarios.  A  clearer  quantitative  analysis  can  be 
achieved  using  the  normal  mode  representation. 

4-1.1  Normal  Mode  Representation 

At  the  receiving  array  we  measure  the  pressure  fields  from  the  target  and  the 
reference  source  as  the  latter  is  deployed  to  various  depths;  while  at  the  virtual 
array  we  compare  the  calculated  pressure  field  from  a  test  source  with  the  holo¬ 
graphic  image.  Therefore,  the  reconstructed  field  depends  mainly  on  the  fidelity 
of  the  holographic  image,  and  the  fidelity  of  the  field  from  the  test  source  as  cal¬ 
culated  along  the  virtual  array. 

We  assume  that  the  virtual  array  is  very  long  and  adequately  sampled  and 
that  the  receiving  array.  Fig.  4-2,  is  long  enough  to  adequately  measure  most  of 
the  modes,  especially  the  lower  order  modes.  Then  the  holographic  image  is 
accurately  estimated  along  the  full  length  virtual  array,  and  the  field  from  the 
test  source  is  adequately  sampled  by  the  long  virtual  array.  As  a  result  we  expect 
no  major  degradation  in  the  reconstructed  image.  Thus,  we  will  be  able  to  recon¬ 
struct  the  wavefront  from  the  limited  information  in  the  original  hologram 
(receiving  array)  almost  perfectly,  because  the  full  length  virtual  array  cancels 
out  the  media  effects  for  all  the  measured  paths,  as  shown  in  Fig.  4-3.  As  the 
receiving  array  becomes  very  short,  the  lower  order  modes  are  poorly  resolved, 
adversely  affecting  the  resolution  of  the  reconstructed  image. 
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FIGURE  4-2.  Diagram  of  the  different  modes  as  measured  by  a  truncated  receiving  array. 


Shortening  the  virtual  array  will  result  in  a  limited  cancellation  of  the 
medium  effects  for  the  measured  paths.  In  Fig.  4-4,  paths  1  and  2  from  the  target 
are  completely  phase  compensated  by  corresponding  paths  from  the  reference 
array.  Paths  3  and  4  are  partially  phase  compensated  by  corresponding  paths 
from  the  reference  array,  path  3  propagates  over  a  range  Xj  from  the  reference 
array  with  no  phase  compensation,  and  path  4  propagates  over  range  x2  with  no 
phase  compensation.  Finally  path  5  is  not  phase  compensated.  Since  the  effect  of 
the  noncompensated,  or  partially  compensated  paths,  can  be  quite  large  for  long 
range  propagation,  virtual  array  truncation  may  result  in  a  holographic  image 
with  distorted  phase.  Also,  the  shorter  virtual  array  will  poorly  resolve  the  lower 
order  modes  which  carry  quite  a  bit  of  information.  Thus  the  combination  of  not 
phase  compensating  certain  paths  and  poorly  resolving  certain  modes  will  result 
in  a  severe  degradation  of  the  reconstructed  image  even  if  the  receiving  array 
spans  the  whole  water  column. 


FIGURE  4-3.  Schematic  representation  of  the  HAP  algorithm  with  a  truncated  receiving 
array.  All  the  paths  from  the  target  are  calibrated  by  the  virtual  array. _ 
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FIGURE  4-4.  Schematic  representation  of  the  FLAP  algorithm  with  a  truncated  virtual 
array.  Some  paths  from  the  target  are  not  completely  calibrated  by  the  virtual  array. 
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Normal  mode  theory  yields  additional  understanding  of  the  different  roles 
that  the  truncated  receiving  and  virtual  arrays  play  in  degrading  the  recon¬ 
structed  image.  Comparing  the  measured  and  estimated  (calculated)  modes  asso¬ 
ciated  with  both  the  receiving  and  virtual  arrays  will  lead  to  a  quantitative  and 
qualitative  evaluation  of  the  degradation  due  to  the  different  schemes  of  array 
truncation. 

Eq.  (2-32)  shows  the  pressure  at  range  r  and  depth  z  from  a  point  source  at 
range  r,  and  depth  ztin  a  range  dependent  adiabatic  medium,  and  is  repeated 
here  for  convenience: 


o  1/2  • 

2  ,  -in/ 4 


P(,r,z-.rs,zs)  =  (-)  e‘  Y.Vm(r,  z)Vm*  (rs,  zs)  Um(r)] 


-1/2 


m 


x  exp 


L  0 


(4-9) 


where  \|tm(r,z)  is  the  eigenfunction  of  the  mth  mode  at  range  r  and  depth  z,  and  £m 
is  the  horizontal  eigennumber  of  the  mth  mode.  Substituting  Eq.  (4-9)  into  Eq. 
(4-3)  we  get  the  modal  representation  of  the  depth  dependent  holographic  image, 


H(rR,zR)  =  H(rR,zR) 


=£1 

R  m  n 


Vm*(rt,  zt)Vn(rR,  zR) 


8 


exp 


x 


r,  rR 

-i^mCr)d-r-il^Cr,d-r 
L  0  0 


J (r0’z>*dz  * 


(4-10) 


When  h  =  h  the  second  term,  Eq.  (4-10),  will  be  zero  and  H  (zR)  =  H  ( zR ) .  The 
other  scenario  is  when  §<h<h\  then  the  second  term,  Eq.  (4-10),  represents  a 
degradation  factor  subtracted  from  the  original  (nonperturbed)  holographic 
image. 


The  normal  mode  representation  of  the  reconstructed  image,  Eq.  (4-6),  is 
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exp 


Q(r’.z')  = 

m  n  l 


-*  |  J “  J  ^  <r' "  r*> 


xVm(rr  2t)  (r',z') 


8nZJvm  (r0>  z>  Ofc  *)  dz  +  6mn  J  Vjl*  (/>  lj?)  (/>  *jj)  dz* 

h  hR 

h  hR 

-Jvm  (r0’  z)  vn*  (''O’  z)  rfz  J  Vn*  <r*’  ZR}  %  {rR>  ZR)  dzR 


hR 


(4-11) 


,  4  (2/7c)3/2  -in/4 

where  A  =  .  . -^=e 


hR 

bnl  =  1  %.*  (r/?’  Z*}  V/  (r/?’  z/?>  ^  *  <4'12) 

0 

and 

h 

bmn  =  K,  (r0>  Z)  (r0’  z>  Jz  *  (4'12)a 

0 

5  is  the  Kronecker  delta  function,  and,  if  we  assume  that  the  measured  and 

mn 

calculated  (estimated)  eigenfunctions  along  the  virtual  array  are  the  same,  then 
Eq.  (4-12)  approximates  5n/.  This  representation,  Eq.  (4-11),  leads  us  to  analyze 
three  cases: 
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L  No  Truncation 

Both  the  receiving  and  reference  array  span  the  ocean  column  from  surface 
to  bottom,  h  =  h  and  Hr  =  hR\  then  the  three  terms  in  brackets,  Eq.  (4-11),  will 
be  zero,  and  no  truncation  effect  is  in  place. 

n.  Truncated  Receiving  Array 

The  receiving  array  stretches  from  the  surface  of  the  ocean  to  a  certain  depth 
while  the  reference  array  spans  the  water  column  from  surface  to  bottom,  h<h 
and  hR  =  hR;  see  Fig.  4-3.  Then  both  the  second  and  third  terms  in  brackets,  Eq. 
(4-11),  are  equal  to  zero,  and  assuming  that  ( rR ) ,  the  reconstructed 

image  is 


exp 


Q  (rf,  z') 


e(r',n-A£X- 

m  l 


-M 


r,  rR 

J (?)  dr  ~  j  S/  d~r  ~  4/  "  r/?) 
0  0 


xvm(r/5 (r',z')\ym(rQ>z)\ \f*  ( rQ,z)dz  . 

h  (4-13) 


For  an  appropriately  sampled  long  receiving  array  we  may  have  an  almost 
complete  set  of  measured  modes,  and  as  a  result  the  eigenfunctions  under  the 
integral  in  Eq.  (4-13)  are  quasi-orthonormal;  this  means  that 


h  r1+E 

JVOT(z)V/*  (z)dz=  |  em/ 
0  ^  0 


for  m  -  l 

for  m  =  l±L 
else 


(4-14) 


where  e  is  a  small  number,  and  Lisa  series  of  integers,  L  =  { 1, 2,  3,  ... }  .  The 
length  of  L  depends  on  the  size  of  the  array,  as  illustrated  in  Fig.  4-5,  which 
shows  the  correlation  of  the  10th  mode  with  the  other  modes  for  various  array 
lengths. 
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For  a  sufficiently  long  array,  only  the  diagonal,  sub-diagonal,  and  super¬ 
diagonal  elements  of  the  eigenfunctions  correlation  matrix  have  significant  val¬ 
ues.  Therefore,  the  residual  integral  in  Eq.  (4-13)  is 

h 

jVm(z)V/*  (z)dz  = 

h 

where  we  assume  that  the  eigenfunctions  are  orthonormal  over  the  full  water 
column.  Thus,  the  residual  integral  in  Eq.  (4-13)  is  very  small.  The  phase  of  the 
error  term,  at  ranges  comparable  with  the  target  range,  is 


,  e 

Zml 

0 


for  m  =  l 
for  m  =  l±  1 
else 


(4-15) 


<Pm,  =  i^iCr)^-ii(r'-rR) 


Assuming  that  rt>rR  then 


rR  rt  rR 

<Pm,  =  0)^-1, 

o  rR  0 

rR  ri 

=  J« m <?)- 5, <?)]^+j5m (?><*- 4/ ('■'-'*>  • 

0  rR 

Introducing  the  mean  modal  value 


(4-16) 


(4-17) 


(4-18) 


Eq.  (4-17)  can  be  written  as 

rR 

<Pm,  =  |  - %, Cr)]<Tr  +  ^„ (r-rR)  . 

0 


(4-19) 


Scanning  in  the  neighborhood  of  the  target,  where  r'  =  r(±  Ar,  and  A r  is  very 
small  range,  then  the  phase  of  the  error  term  is 
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rR 

0 

rR 

=  |  I5„(?) -$,m<rr+  (£-$,)  (r,-rfi)+|,Ar  . 

0  (4-20) 

In  the  limit,  as  Ar  shrinks  to  zero,  the  phase  of  the  error  term  becomes 


rR 

<Pmr  1  C-y  (r,-r,)  .  (4-21) 

0 

From  the  above,  we  conclude  that  the  coupling  of  the  neighboring  modes 
contributes  to  the  error  term  in  Eq.  (4-13)  with  a  scaled  modal  amplitude  Eml 
and  phase  cp^.  Thus,  the  error  term  is 


",{P  m, 


E  (r\  z')  =  (r,,  2,)  V,  (>•',  z')  e„(  .  (4-22) 

m  '  Kv, 


For  widely  separated  modes,  the  mode  coupling  integral  is  very  small,  em/  =  0, 
see  Fig.  4-5;  therefore,  the  error  term  is  negligible,  E  ~  0.  For  the  neighboring 
modes,  the  difference  between  the  modal  horizontal  wavenumbers,  is 

very  small,  on  the  order1  of  10-4  to  10-6  (the  closer  the  modes  are,  the  smaller 
the  modal  difference  is);  see  Fig.  3-4.  Assuming  that  =  ^  between  r(  and  rR, 
then  the  main  contribution  to  the  phase  fluctuation,  when  r'  ~  r{,  is  from  the  first 
term  in  Eq.  (4-21).  The  modal  difference  in  the  first  term  is  integrated  along  the 
receiver/reference  separation,  which  is  on  the  order  of  hundreds  of  kilometers, 
while  in  the  second  term,  the  modal  difference  is  integrated  over  the  target/ref¬ 
erence  separation,  which  is  on  the  order2  of  tens  of  kilometers.  Therefore,  a 
small  change  in  the  difference  in  horizontal  wavenumbers  produces  a  significant 


1.  These  values  are  taken  from  the  adiabatic  ocean  model  that  is  presented  in  chapter  3. 
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fluctuation  in  the  phase,  and  the  error  terms  become  oscillatory,  resulting  in  a 
noncoherent  summation  over  all  modes.  Therefore,  the  total  error  term  tends  to 
zero  and  the  image  degradation  effect  will  be  reduced. 


An  excessively  short  receiving  array  will  result  in  poorly  resolved  lower 
order  modes  which  propagate  close  to  the  sound  speed  channel.  Fig.  4-2.  Only 
the  higher  order  modes  will  be  adequately  sampled  and  well  represented  by  the 
measurements.  Also,  for  a  receiving  array  much  shorter  than  the  depth  of  the 
channel  axis,  the  quasi-orthonormal  property  of  the  eigenfunctions  will  not  be 
satisfied  over  the  receiving  array,  but  it  will  be  sustained  over  the  complemen¬ 
tary  depth  of  the  receiving  array,  h  to  h. 


h  #  r1+e 

Jvw(*)ty*  (z)dz  =  \  eml 

h  1  0 

Thus  the  error  term  in  Eq.  (4-13)  will  be 


for  m  =  l 

for  m  =  l±L 
else 


(4-23) 


EIK(r',z')=A  XXT - - 


=' (v  2,)  vf  *')  (Sm,  +  eMj)  (4-24) 


where  is  the  Kronecker  delta.  Considering  Eq.  (4-17),  the  error  term 
becomes 


2.  For  the  ocean  model  considered  in  the  previous  chapter,  rR  -  500  km  and  rt  -  rR  -  20  km .  Assuming  an 
average  wavenumber  difference  on  the  order  of  10-5 ,  then,  the  first  term  in  Eq.  (4-21)  is  on  the  order  of  5 
rad,  and  the  second  term  is  on  the  order  of  0.2  rad. 
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-1/2 


ETR  (r'f  z')  =  V  Vm*  (r'«  z'>  [$m  (r/)  <r*>  5  J  * 

m 


exp 


■r, 

-/• 

L  ^jj 


<r'-rs) 


exp  [— i cp  ,]  * 

+  /1I£  I  "  ■  <-r'-2">tml 

m  ‘  Jtm(r,)i,(rK)i, 


(4-25) 


Because  the  first  term  is  equal  to  Q  ( r\  z')  we  can  write 


E™  («•',  2')  =  Q  (r',  2')  +  4  TYw (<>  z,>  V*  (' *')  x 

m  l 

[A»<rl>Mr*>$,]  e„,,  • 

Substituting  Eq.  (4-26)  in  Eq.  (4-13)  we  get 


(4-26) 


★  r~  *  -1  1/2 

e  (/•',  2-)  =  z,>  v/  (<■'. z'>  [tM  (<•,)  $,<'*>  x 

m  / 

«p[-i9m,]emJ  .  (4  2?) 

Therefore,  when  the  receiving  array  is  excessively  short,  the  quasi-orthogonality 
property  does  not  hold,  and  the  reconstructed  image  has  large  error;  therefore, 
the  algorithm  will  fail. 

HI.  Truncated  Virtual  Array 

The  receiving  array  spans  the  water  column  of  the  ocean  from  surface  to 
bottom,  h  =  h,  while  the  reference  source  is  deployed  from  surface  to  a  certain 
depth  in  the  water  creating  a  limited  length  virtual  array,  hp  <  hR  (see  Fig.  4-4). 
The  first  and  third  terms  in  brackets,  Eq.  (4-11),  are  equal  to  zero  and  the  recon¬ 
structed  image  is 
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exp 


Q(r',z')  = 


Q(r',z')-AZX- 

m  l 


-i\%m  (r)  dr—i  j  i,n  (?)  dr  -  (/ -  rR) 

.00 


Vm  <rr  z,)  V*  <r\  z')  j  V„*  (zR,  zR)  V,  (r„,  zfi)  dzR  . 

h  (4-28) 


A  similar  argument  as  in  the  previous  case,  regarding  the  quasi-orthonormal 
property  of  the  eigenfunctions  over  a  truncated  array,  is  valid.  Comparing  Eq. 
(4-28)  with  Eq.  (4-13)  and  assuming  that  the  water  column  at  the  virtual  array  is 
of  the  same  depth  as  that  at  the  receiving  array,  hR  =  h,  and  the  truncated  virtual 
array  is  of  the  same  length  as  the  truncated  receiving  array,  hR  =  h,  then  Eq.  (4- 
28)  is  similar  to  Eq.  (4-13),  except  for  the  phase  term  and  the  eigenfunction  cor¬ 
relation  matrix.  In  Eq.  (4-13)  the  eigenfunction  correlation  matrix  is  found  for 
the  measured  eigenfunctions  along  the  receiving  array;  however,  in  Eq.  (4-28) 
the  cross  correlation  matrix  between  the  measured  and  the  calculated  (esti¬ 
mated)  eigenfunctions  is  found.  The  measured  eigenfunctions,  \|/m  ( rR ,  zR ) ,  can 
be  found  from  the  modal  decomposition  of  the  measured  pressure  field,  from  the 
reference  source  at  (rR,zR),  along  the  receiving  array.  The  estimated  eigenfunc¬ 
tions,  (rR,  zR ) ,  are  found  using  the  assumed  sound  speed  profile  at  the  refer¬ 
ence  location. 


Assuming  an  agreement  between  the  measured  and  estimated  eigenfunc¬ 
tions,  \\fm  (rR,  zR)  =  \j fm  (rR,  zR) ,  the  correlation  matrices  in  Eqs.  (4-13)  and  (4- 
28)  are  the  same  and  the  algorithm  performance  with  depth  in  both  cases  will  be 
equivalent.  Otherwise,  a  poor  estimate  of  the  eigenfunctions  along  the  virtual 
array  will  result  in  performance  degradation  in  the  vertical  (defocusing  with 
depth). 

The  difference  in  the  phase  term  for  the  two  truncation  scenarios  is  an  indi¬ 
cation  of  different  algorithm  behavior  (focusing)  in  range.  The  phase  of  the 
error  term  in  Eq.  (4-28)  is 
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0  0 

assuming  that  rt>rR\  then  the  phase  of  the  error  term  can  be  written  as 


(4-29) 


KirXfr+l^rXTr 


-\t>mCr)(fr-\l{r'-rR) 


(4-30) 


Scanning  in  the  neighborhood  of  the  target,  where  r'  =  r,±Ar,  and  A r  is  very 
small  range,  and  using  the  modal  mean  value,  Eq.  (4-18),  then  Eq.  (4-30)  can  be 
written  as 


<P ml  =  ^rrrR^-^rt  +  Ar-rR^  •  <4'31> 

In  the  limit  as  A r  shrinks  to  zero,  the  phase  of  the  error  term  becomes 

<Pm,  =  •  (4-32) 

From  the  above,  we  conclude  that  the  coupling  of  the  neighboring  modes 
contributes  to  the  error  term  in  Eq.  (4-28)  with  a  scaled  modal  amplitude  zml 
and  phase  <pm/.  Therefore,  the  error  term  can  be  written  as 


ETV(.r',z')  (z')£m/  •  <4'33> 

m '  Kifii 

Notice  that  the  phase  of  the  error  term  is  a  function  of  the  target/reference 
separation  only.  For  the  far  apart  modes,  the  scaled  modal  amplitude  is  very 
small,  e  j  =  0,  and  as  a  result  the  error  term  is  negligible,  E  =0.  For  neighbor¬ 
ing  modes,  and  assuming  =  then  is  very  small,  on  the  order  of 

10~4  to  10~6.  Therefore,  small  changes  in  the  horizontal  wavenumber  difference, 
1^— will  not  result  in  drastic  fluctuations  in  the  phase,  since  these  changes 
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are  integrated  over  the  target/reference  separation  which  is  on  the  order  of  10 
km.  This  will  lead  to  a  coherent  summation  of  most  of  the  terms  in  Eq.  (4-33), 
and  as  a  consequence  the  total  error  can  be  significant,  leading  to  severe  degra¬ 
dation  (defocusing)  of  the  reconstructed  image. 

Normalizing  the  magnitude  of  the  truncation  error  term  by  the  magnitude  of 
the  truncation-free  reconstructed  field,  we  get 


ETR(r',  z') 

\Q(r',z')\ 


for  the  truncated  receiving  array,  and 


(4-34) 


ETV(r',  z') 
\Q(r’,z')\ 


(4-35) 


TR 

for  the  truncated  virtual  array.  E  is  the  error  due  to  truncating  the  receiving 


array,  Eq.  (4-22),  and  E 
33). 


TV 


is  the  error  due  to  truncating  the  virtual  array,  Eq.  (4- 


Fig.  4-6  shows  the  normalized  magnitude  of  the  error  terms  as  a  function  of 
range,  r'  —  r  .  The  truncated  receiving  array  and  the  truncated  virtual  array  are  of 
the  same  length  (1500  m);  see  Fig.  4-5  for  the  residual  integrals  in  the  error 
terms  for  this  array  length.  The  ocean  model  is  the  same  as  that  used  in  chapter 
3  with  exponentially  stratified  SSP.  For  all  ranges  considered  (0  to  5000  meters 
from  the  target),  we  notice  that  the  truncated  virtual  array  has  a  larger  error  term 
(solid  line)  than  a  comparable  truncated  receiving  array  (doted  line).  Also,  the 
truncated  virtual  array  is  more  susceptible  to  changes  in  A r,  see  Eq.  (4-31), 
which  is  evident  from  the  oscillatory  behavior  of  the  solid  line. 

It  can  be  concluded  that  the  length  of  the  virtual  array  is  more  important 
than  the  length  of  the  receiving  array  to  the  algorithm’s  performance  and  the 
quality  of  the  reconstructed  image. 
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range  (m) 


FIGURE  4-6.  The  magnitude  of  the  normalized  error  term  for  truncated  receiving  array 
(Eq.  4-22)  (doted  line),  and  that  for  truncated  virtual  array  (Eq.  4-33)  (solid  line).  The 
horizontal  scale  shows  the  range  of  the  test  source  from  the  target. 


4-2  Simulation  Results 

The  effect  of  truncating  the  receiving  array,  virtual  array,  or  both  on  the  per¬ 
formance  of  the  HAP  algorithm  is  evaluated  using  a  high  order  PE  code,  FEPE, 
in  a  range  dependent  ocean  with  exponentially  stratified  sound  speed  profile, 
Munk  profile,  where  the  axis  depth  is  1000  meters  and  the  range  dependency  is 
modeled  as  a  linear  variation  of  the  axial  sound  speed  with  range.  The  slope  of 
this  variation  is  taken  to  be  one  meter  per  second  per  100  km. 

Under  ideal  conditions,  full  length  receiving  array  and  full  length  virtual 
array,  we  showed  that  the  holographic  array  processor  produces  a  tight  focus  at 
the  target  location  with  relatively  low  level  side  lobes.  Fig.  3-15. 

4-2.1  Truncated  Receiving  Array  Simulation 

We  assume  that  the  length  of  the  receiving  array  is  limited.  It  extends  from 
the  surface  of  the  ocean  to  a  certain  depth  less  than  the  total  depth  of  the  water 
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column;  at  the  same  time,  the  reference  source  spans  the  whole  water  column, 
creating  a  full  length  virtual  array,  as  illustrated  in  Fig.  4-7.  The  target  is  located 
at  25  km  from  the  virtual  array  which  in  turn  is  located  at  450  km  from  the 
receiving  array. 


r 


FIGURE  4-7.  Schematic  diagram  of  the  ocean  model  with  range  dependent  Munk  SSP 
and  truncated  receiving  array.  _ 


Fig.  (4-9)  shows  color  images  of  the  magnitude  of  the  reconstructed  fields, 
Q  ( r z) ,  for  various  receiving  array  lengths.  The  effect  of  truncating  the 
receiving  array  on  the  performance  of  the  HAP  algorithm  is  evident.  Fig.  4-10 
shows  horizontal  slices,  at  the  target’s  depth,  of  the  normalized  power  of  the 
reconstructed  images,  Q  (/,  zt) ,  for  various  receiving  array  lengths.  It  gives  a 
clear  picture  of  the  effect  of  truncating  the  receiving  array  on  the  width  of  the 
main  lobe  (focus)  and  the  level  of  the  side  lobes  (ambiguity).  We  see  that  when 
the  receiving  array  covers  the  upper  half  of  the  water  column,  2500  meters,  the 
algorithm  does  not  suffer  from  any  major  degradation;  it  produces  a  sharp  focus 
on  the  target  with  negligible  rise  in  the  side  lobes’  levels.  As  we  reduce  the 
length  of  the  receiving  array  down  to  a  depth  comparable  to  the  depth  of  the 
channel  axis,  1000  meters  for  this  ocean  model,  the  major  effect  on  the  algo- 
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rithm  performance  is  a  relatively  slight  increase  in  the  side  lobes,  but  the  focus 
is  still  sharp  at  the  target.  When  we  implement  a  receiving  array  that  extends 
from  the  surface  to  a  depth  shallower  than  the  depth  of  the  channel  axis,  the 
algorithm  fails  to  produce  a  sharp  focus  and  the  side  lobes  rise  to  high  levels, 
comparable  to  the  main  lobe’s  level. 

4-2.2  Truncated  Virtual  Array  Simulation 

The  receiving  array  spans  the  water  column  of  the  ocean.  The  reference 
source  is  deployed  to  a  certain  depth,  less  than  the  total  depth  of  water  column 
at  the  reference  site,  creating  a  truncated  virtual  array.  The  ocean  model  and  the 
separations  between  the  target,  virtual  array,  and  receiving  array  are  the  same  as 
in  the  previous  section;  see  Fig.  4-8. 


FIGURE  4-8.  Schematic  diagram  of  the  ocean  model  with  Munk  SSP  and  truncated 
virtual  array. 


Fig.  4-11  shows  color  images  of  the  magnitude  of  the  reconstructed  fields, 
Q  (r\  zt ) ,  for  various  virtual  array  lengths.  The  effect  of  the  virtual  array  length 
on  the  performance  of  the  HAP  algorithm  is  clear.  Fig.  4-12  shows  horizontal 
slices,  at  the  target’s  depth,  of  the  normalized  power  of  the  reconstructed 
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images,  Q  (r',  z() ,  for  various  virtual  array  lengths.  As  in  Fig.  4-10,  this  plot 
gives  a  clear  picture  of  the  effect  of  truncating  the  virtual  array  on  the  width  of 
the  main  lobe  (focus)  and  the  level  of  the  side  lobes  (ambiguity).  Deploying  the 
reference  source  down  to  half  the  water  column  will  create  a  virtual  array  of 
length  2500  meters;  in  this  case  the  algorithm  does  not  suffer  from  any  major 
degradation,  and  a  sharp  focus  at  the  target  is  produced  with  low  level  side 
lobes.  As  we  deploy  the  reference  source  to  shallower  depths,  we  notice  widen¬ 
ing  in  the  main  lobe  (defocusing)  and  a  rise  in  the  side  lobes  (ambiguity).  This 
degradation  gets  worse  and  worse  with  decreasing  size  of  the  virtual  array,  down 
to  a  point  where  the  size  of  the  virtual  array  is  comparable  to  the  depth  of  the 
channel  axis;  then  the  main  lobe  becomes  very  fat  (smeared  focus)  and  the  side 
lobes  rise  to  or  exceed  the  level  of  the  main  lobe  (strong  ambiguity). 

4-2.3  Conclusion 

The  results  of  the  numerical  simulation  agree  with  the  theoretical  findings 
which  were  presented  in  the  previous  section.  As  a  measure  of  comparison 
between  the  effect  of  truncating  the  receiving  array  and  the  virtual  array  on  the 
performance  of  the  algorithm,  we  plot  the  range  band  width  of  the  half  power 
level  of  the  main  lobe  as  a  function  of  the  array  length.  Fig.  4-13;  the  array 
length  is  measured  in  terms  of  the  acoustic  wave  length.  It  is  evident  from  this 
plot  that  the  length  of  the  virtual  array  is  crucial  to  the  performance  of  the  algo¬ 
rithm  and,  therefore,  to  the  quality  of  the  reconstructed  image. 

The  simulation  was  done  for  an  acoustic  frequency  of  50  Hz  (acoustic  wave 
length  30  meters).  We  see  that  as  the  length  of  the  virtual  array  becomes  close  to 
the  depth  of  the  channel  axis  (1000  meters),  the  width  of  the  main  lobe  becomes 
a  significant  portion  of  the  target/  reference  separation,  rendering  the  algorithm 
ineffective.  For  example,  when  the  virtual  array  extends  from  the  surface  to  a 
depth  of  1500  meters,  the  half  power  range  band  width  is  about  4  km;  this  is 
equivalent  to  16%  of  the  target/reference  separation.  If  we  deploy  the  reference 
to  a  yet  shallower  depth,  such  as  1250  meters,  then  the  half  power  range  band¬ 
width  is  about  9  km  which  is  equivalent  to  36%  of  the  target/reference  separa¬ 
tion.  This  drastic  variation  in  the  algorithm  performance  is  not  present  with  a 
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change  in  length  of  the  receiving  array,  where  the  half  power  range  band-width 
is  less  than  1  km  which  constitutes  about  4%  of  the  target/reference  separation. 

It  is  noticeable  from  the  plots  that  the  level  of  the  side  lobes  is  more  suscep¬ 
tible  to  the  length  of  the  virtual  array  than  to  the  length  of  the  receiving  array. 
Therefore,  we  conclude  that  the  length  of  the  virtual  array  is  more  crucial  to  the 
focusing  and  ambiguity  performance  of  the  HAP  algorithm. 


depth  (km) 
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FIGURE  4-9.  Magnitude  of  Q  as  a  function  of  range  and  depth  where  the  receiving  array 
length  varies:  (a)  1000  m,  (b)  1250  m,  (c)  1500  m,  (d)  1750  m,  (e)  2000  m,  (f)  2500  m. 
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FIGURE  4-11.  Magnitude  of  Q  as  a  function  of  range  and  depth  where  the  virtual  array 
length  varies:  (a)  1000  m,  (b)  1250  m,  (c)  1500  m,  (d)  1750  m,  (e)  2000  m,  (f)  2500  m. 


Ar  (m) 


FIGURE  4-13.  The  half  power  range  band- width  as  a  function  of  the  array  length. 


CHAPTER  5 


Ambient  Noise 

The  performance  of  the  HAP  algorithm  in  the  presence  of  a  noise  field  can 
be  analyzed  through  the  knowledge  of  the  time-space  statistics  of  the  field 
across  the  receiving  array. 

The  background  acoustic  noise,  or  ambient  noise,  establishes  a  lower  limit 
on  the  intensity  of  useful  acoustic  signals.  Many  useful  signals  in  the  ocean  are 
noise-like  in  character.  For  instance,  the  presence  of  many  ships,  randomly  dis¬ 
tributed  over  the  ocean  surface,  results  in  a  component  of  ambient  noise 
ascribed  to  “distance  shipping”  or  “ship  traffic”.  However,  the  noise  produced 
by  a  single  nearby  ship  is  easily  identified  and  localized  and  is  therefore  treated 
as  an  acoustic  signal  rather  than  as  part  of  ambient  noise. 

5-1  Sources  of  Ambient  Noise 

Ambient  noise  can  be  divided  into  three  categories;  natural,  biological,  and 
human-made  sources.  Each  noise  source.  Fig.  5-1,  exhibits  different  directional 
(spatial)  and  spectral  (temporal)  characteristics.  Among  the  natural  sources  of 
noise  are  seismic  disturbances,  agitation  of  the  sea  surface  by  wind,  and  thermal 
activity  of  the  water  molecules.  Biological  sources  include  whales,  snapping 
shrimp,  croakers,  dolphins,  and  various  other  fish  and  ocean  mammals.  The 
principal  human-made  component  of  ambient  noise  is  the  sound  generated  by 
distant  (over  hundreds  of  nautical  miles)  shipping,  and  off-shore  drilling. 

It  is  evident  that  ambient  noise  levels  are  highly  dependent  on  geographic 
location,  acoustic  transmission  characteristics,  season  of  the  year,  and  weather. 
Even  political  events,  such  as  the  closing  or  opening  of  a  major  canal  or  port  and 
opening  or  closing  of  trade  and  tourism  between  different  countries,  have  a  mea¬ 
surable  effect  on  the  ambient  noise  through  their  influence  on  global  shipping 
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Table  (5-1)  shows  the  dominant  noise  sources  in  the  different  spectrum 
bands.  The  very  low  frequency  ambient  noise  (f  <  10  Hz)  is  generally  not  of 
great  interest  in  the  design  of  underwater  acoustic  systems. 
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FIGURE  5-1.  Conceptual  diagram  showing  some  of  the  sources  of  deep-water  ambient 
noise  as  measured  at  a  hydrophone.  _ 


Table  5-1:  The  dominant  noise  sources  in  the  different  spectrum  bands 


frequency 

dominant  sources 

comment 

f  <  10  Hz 

seismic  and  ocean  turbu¬ 
lent  pressure  fluctuations 

Very  low  frequency. 

10<f<  100  Hz 

distant  shipping  and  off¬ 
shore  drilling 

Low  frequency.  Broad  peak 
around  30  Hz. 

100  <f<  1000  Hz 

wind 

Mid  frequency. 

Broad  peak  around  500  Hz 

f  >  50  kHz 

thermal  agitation  of  the 
water  molecules 

High  frequency.  Increases  at  a 
rate  of  6  dB/octave  as  fre¬ 
quency  increases 

The  low-frequency  ambient  noise  arriving  from  great  distances  depends  pri¬ 
marily  on  refraction  in  the  ocean  sound  channel.  Paths  requiring  multiple  bottom 
reflections  endure  large  attenuation  and  therefore  do  not  propagate  very  far. 
Hence,  high  ambient  noise  levels  in  the  low  frequency  region  generally  occur  in 
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regions  where  water  depths  and  sound  speed  profiles  support  convergence  zone 
transmission.  For  typical  sound  speed  profiles  in  the  deep  ocean,  the  vertical 
angle  of  arrival  at  a  receiving  hydrophone,  for  acoustic  rays  that  do  not  intersect 
either  the  ocean  surface  or  bottom,  is  limited  approximately  to  ±15°  relative  to 
the  horizontal.  Fig.  5-2  shows  a  typical  angular  noise  level  distribution  in  the 
vertical  plane  for  the  low-frequency  band.  The  intensity  level  is  in  dB  with 
respect  to  one  pPa  in  a  one  Hz  band  per  steradian  solid  angle. 


FIGURE  5-2.  A  typical  angular  noise  level  distribution  in  the  vertical  plane. 


In  the  midfrequency  band,  ambient  noise  is  dominated  by  the  effects  of  wind 
activity  on  the  sea  surface.  In  very  approximate  terms  [20],  the  ambient  noise 
level  in  this  band  increases  about  5  dB  for  every  doubling  of  wind  speed  with  a 
wide  spectral  peak  around  500  Hz  and  decreases  at  a  rate  of  -5  to  -6  dB  per 
octave  at  higher  frequencies.  Because  of  the  propagation  and  geometric  proper¬ 
ties  of  the  surface  generated  ambient  noise,  the  intensity  level  of  this  noise  tends 
to  be  higher  in  the  vertical  direction  than  in  the  horizontal  direction. 

Average  representative  ambient  noise  spectra  for  different  conditions  are 
adequate  for  qualitative  characterization.  Fig.  5-3  shows  the  Wenz  curves  [21] 
for  different  conditions  of  shipping  and  wind  speed. 
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When  more  than  one  source  of  noise  is  present,  the  effective  noise  back¬ 
ground  is  obtained  by  summing  the  intensities  of  the  contributing  sources.  When 
using  noise  levels  (specified  in  units  of  spectrum  level,  dB  relative  to  1  pPa), 
this  summation  process  is  easily  accomplished  using  the  ‘power  summation’ 
operator,  denoted  by  the  symbol  ©,  and  defined  as 


©  =  lOlog  10L,/1°  (5-1) 

i  =  l 

where  Li  is  the  level  of  i-th  noise  source  in  dB  and  n  is  the  number  of  contribut- 

' _ _  _  ! _ _  rT'*l_  !  -  oti  tialii  O  *■  A  «  CA  IaITaI  /  \  t  A 


For  example,  in  Fig.  5-3  at  frequency  100  Hz  under  conditions  of  moderate 
shipping  and  sea  state  6,  one  would  obtain  noise  levels  of  about  69  dB  and  71 
dB,  respectively.  The  noise  level  is  obtained  from  Eq.  (5-1)  as 


69dB&71dB  =  lOlog  [  1069  +  1071]  =  73 dB 


(5-2) 
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which  is  2  dB  higher  than  the  noise  level  due  to  surface  weather  alone. 

5-2  Depth  Dependence 

The  ambient  noise  intensity  is  a  very  weak  function  of  the  receiver  depth 
when  operating  in  the  low  frequency  band,  f  <  500  Hz.  For  higher  frequencies 
the  depth  dependence  becomes  significant  due  to  absorption. 

Measurements  of  the  depth  dependence  of  low-frequency  ambient  noise 
were  made  by  Morris  (1978)  [22]  in  the  Northeastern  Pacific  Ocean.  Hydro¬ 
phones  were  suspended  from  the  research  platform  FLIP  (Floating  Instrument 
Platform).  Fig.  5-4  shows  the  sound  speed  profile,  hydrophone  depths  and  aver¬ 
age  noise  profiles  in  one-third-octave  bands  for  this  experiment.  There  is  a 
decrease  of  noise  level  with  depth  at  low  frequencies,  with  a  smaller  decrease 
with  depth  at  500  Hz  as  wind  noise  overcomes  the  dominance  of  shipping  noise. 
Below  the  critical  depth,  where  the  sound  speed  equals  that  at  the  surface,  the 
fall-off  with  depth  is  steeper  as  the  bottom  is  approached.  This  is  the  result  of 
the  loss  of  refracted  sound  energy.  A  hydrophone  below  the  critical  depth 
receives  energy  only  through  reflection  from  the  surface,  the  bottom,  or  both, 
therefore  the  received  energy  from  a  distant  source,  at  these  depths,  will  be 
highly  attenuated  due  to  interaction  with  the  surface,  the  bottom,  or  both. 
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FIGURE  5-4.  (a)  SSP  and  hydrophone  distribution  in  the  water  column,  (b)  Variation  of 
ambient-noise  level  with  depth  at  the  1978  experiment  in  the  Northeast  Pacific  Ocean. 
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5-3  Directionality 

As  a  first  order  approximation,  the  noise  field  in  the  ocean  might  be  consid¬ 
ered  to  be  isotropic  in  nature,  that  is,  uniform  in  all  directions,  both  horizontal 
and  vertical.  This  approximation  is  not  accurate  yet  it  is  useful  for  qualitative 
analysis. 

Measurements  have  shown  that  the  ambient  noise,  in  general,  is  not  isotro¬ 
pic.  The  vertical  directionality  of  ambient  noise  was  measured  [23]  at  frequen¬ 
cies  of  112  and  1414  Hz  using  a  300  ft.  long  40  element  vertical  array,  located 
south  of  Bermuda.  Fig.  5-5  presents  polar  plots  of  the  ambient  noise  intensity 
per  unit  solid  angle  N( 0)  arriving  at  a  bottomed  hydrophone  as  a  function  of 
vertical  angle  0.  At  112  Hz,  more  noise  appears  to  arrive  at  the  hydrophone  from 
the  horizontal  than  from  the  vertical;  this  difference  diminishes  with  increasing 
wind  speed.  At  1414  Hz  the  opposite  is  true  in  that  more  noise  arrives  from  the 
overhead  than  horizontally;  this  effect  increases  with  increasing  wind  speed. 


112  Hz  1414  Hz 


FIGURE  5-5.  Distribution  of  ambient  noise  intensity  in  the  vertical  plane  at  a  bottomed 
hydrophone  (a)  at  112  Hz,  and  (b)  at  1414  Hz. 
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This  directional  behavior  is  consistent  with  the  view  that  low-frequency 
noise  originates  at  great  distances  and  arrives  at  the  measurement  hydrophone 
via  horizontal  paths,  suffering  little  attenuation,  whereas  high  frequency  noise 
originates  locally  at  the  sea  surface  more  nearly  overhead. 

The  horizontal  or  azimuthal  directionality  of  ambient  noise  can  be  highly 
variable,  particularly  at  low  frequencies.  Shipping  traffic  is  the  dominant  source 
of  noise  at  low  frequencies,  and  the  temporal  and  spatial  variations  in  shipping 
densities  explain  much  of  the  observed  azimuthal  variation. 

5-4  Mathematical  Models 

Mathematical  models  of  noise  in  the  ocean  predict  the  level  and  directional¬ 
ity  (vertical  and  horizontal)  of  noise  as  a  function  of  frequency,  depth,  geo¬ 
graphic  location  and  time  of  year. 

Noise  models  can  be  segregated  into  two  categories;  ambient  noise  models 
and  beam-noise  statistics  models.  Ambient  noise  models  predict  the  mean  levels 
sensed  by  an  acoustical  receiver  when  the  noise  sources  include  surface  weather, 
biology,  and  commercial  activities  such  as  shipping  and  drilling.  Beam-noise 
statistics  models  predict  the  properties  of  low  frequency  shipping  noise  for 
application  to  large-aperture,  narrow  -beam  passive  sonar  systems.  The  latter 
models  use  either  analytic  or  simulation  techniques  to  generate  statistical 
descriptions  of  the  beam  noise.  In  this  context,  beam  noise  is  defined  as  the  con¬ 
volution  of  the  receiver  beam  pattern  with  the  sum  of  the  intensities  from  the 
various  noise  sources.  The  analytical  models  calculate  statistical  properties 
directly  from  the  components  (e.g.  source  level,  propagation  loss),  while  the 
simulation  models  use  Monte  Carlo  techniques. 

Regardless  of  type,  noise  models  consist  of  two  components:  a  transmission 
loss  component,  and  a  noise  level  and  directionality  component.  In  principle,  the 
transmission  loss  can  be  computed  internal  to  the  noise  model  or  it  can  be  input 
externally  from  other  (stand-alone)  model  predictions  or  from  field  measure¬ 


ments. 
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5-4.1  Ambient  Noise  Models 

Ambient  noise  models  consider  variable  densities  of  noise  sources  distrib¬ 
uted  over  large  areas.  This  presents  a  fair  approximation  to  the  generation  of 
wind  noise  as  well  as  distant  shipping  noise.  Consequently,  the  transmission  loss 
calculation  can  be  range  averaged,  as  opposed  to  a  point-to-point  estimate.  This 
greatly  relaxes  the  accuracy  to  which  transmission  loss  must  be  known. 

A  simple  model  of  ambient  noise  in  the  ocean  consists  of  an  infinite  layer  of 
uniform  water  with  a  plane  surface  along  which  the  noise  sources  are  distributed 
and  the  noise  level  is  essentially  independent  of  depth.  A  more  realistic  model, 
however,  would  include  volume  absorption  and  the  effect  of  refraction  and 
boundary  reflections  over  long-range  paths. 

The  low-frequency  component  of  noise  due  to  distant  shipping  can  be  com¬ 
puted  from  knowledge  of  the  density  of  shipping  as  a  function  of  azimuth  and 
range  from  the  receiver  [24],  the  source  level  for  the  radiated  noise  of  each 
generic  type  of  merchant  ship,  and  the  transmission  loss  as  a  function  of  range 
between  the  near  surface  sources  and  the  depth  of  the  receiver.  Then,  contribu¬ 
tions  from  successive  range  rings  centered  about  the  receiver  can  be  summed  to 
obtain  the  level  of  shipping  noise  as  a  function  of  azimuth  at  the  receiver. 

The  high-frequency  component  of  noise  due  to  surface  weather  is  usually 
computed  under  the  assumption  that  it  is  locally  generated  and  isotropic.  Thus, 
only  the  weather  conditions  (sea  state  or  wind  speed)  prevailing  in  the  immedi¬ 
ate  vicinity  need  to  be  considered  in  addition  to  any  localized  rain  showers  or 
biologic  activity. 

Modeling  of  the  vertical  directionality  of  deep-water  ambient  noise  can  be 
approached  by  means  of  a  simple  model  [2].  Consider  a  bottomless,  uniform 
ocean,  without  refraction  or  attenuation,  having  a  surface  covered  with  a  dense, 
uniform  distribution  of  noise  sources.  Furthermore,  let  each  unit  area  of  the  sur¬ 
face  radiate  with  an  intensity  /  (0)  at  one  meter.  Then,  at  point  P,  in  Fig.  5-6  (a), 
the  incremental  intensity  dl  produced  by  a  small  circular  annulus  of  area  dA  at 
horizontal  range  r  is 
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I(Q)dA  I(Q)2%rdr 

dI  =  — J—  =  - J -  *  (5*3) 

But  r  =  Titan  (0) ,  so  that  dr  -  Tisec2  (0)  dQ.  Also,  /  =  Tisec(0).  Substituting 
and  rearranging,  Eq.(5-3)  becomes 


dl  =  2k/ (0)  tan  (Q)dQ  . 


(5-4) 


(a)  (b) 


FIGURE  5-6.  Simple  model  for  the  vertical  directionality  of  ambient  noise;  (a)geometry 
with  straight-line  propagation  paths;  (b)  directional  patterns  for  surface  distribution  of 
monopoles  and  dipoles.  The  dashed  segments  near  the  horizontal  show  the  effect  of 
attenuation  and  refraction  near  the  sea  floor  (A)  and  near  the  deep  sound  channel  axis  (B). 


If  Q  is  a  solid  angle,  then  dQ  =  27tsin  (0)  dQ,  so  that  the  noise  intensity  per  solid 
angle,  N(0),  becomes 


N(Q]  =  =  /  (0)  sec  (0) 


When  I  (0)  =  /Q,  then  Eq.(5-5)  yields 


(5-5) 


N(Q)  =  7osec(0) 


(5-6) 
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which  is  the  ambient  noise  intensity  per  solid  angle  for  a  surface  distribution  of 
monopole  sources.  Using  distribution  of  dipole  sources  (formed  by  the  actual 
source  and  its  image  in  the  sea  surface),  the  intensity  radiated  at  an  angle  0  is 
7(0)  =  /Q cos2  (0) ,  and  the  beam  pattern  of  the  noise  received  at  a  depth  below 
the  surface  is 


IV  (0)  =  /Ocos(0)  .  (5-7) 

At  angles  near  the  horizontal  (0  =  90°)  the  effects  of  attenuation,  refraction  and 
boundary  multipath  prevent  the  intensity,  for  the  monopole  and  the  dipole  mod¬ 
els,  from  going  to  either  infinity  or  zero. 

Using  this  simple  model,  a  receiver  located  within  the  deep  sound  channel 
would  not  be  expected  to  receive  noise  arriving  along  paths  near  the  horizontal 
since  range-independent  ray  tracing  would  show  that  such  ray  paths  do  not  exist. 
The  beam  pattern  of  the  noise  from  monopole  sources  would  therefore  have  a 
maximum  at  an  angle  above  and  below  the  horizontal  (typically  ±10°  to  ±  15°). 
This  conceptual  picture  is  not  always  valid  in  realistic  environments;  specifi¬ 
cally,  more  sophisticated  models  would  consider  the  effects  of  range-dependent 
refraction,  bottom  reflection  and  multipath  arrivals. 

For  low-frequency  noise  sources  (/<  500  77z),  the  agreement  between  avail¬ 
able  observations.  Fig.  5-5,  and  the  findings  of  this  simple  model,  Fig.  5-6,  fur¬ 
ther  suggests  that  a  distribution  of  monopoles  (/  (0)  =  70)  adequately  models 
the  distant  shipping  noise. 

5-5  The  Correlation  Function  for  Ambient  Noise 

The  performance  of  an  array  of  sensors  in  a  noise  field  depends  on  the  cross- 
correlation  function  or  equivalently  the  cross- spectral  density  function  which 
describes  the  second-order  statistical  relationships  between  all  pairs  of  sensors 
in  the  array.  This  function  in  turn  depends  on  the  directional  and  temporal  prop¬ 
erties  of  the  noise  field  and  the  geometry  of  the  array.  A  common  model  for  the 
ambient  noise  assumes  the  noise  to  be  generated  from  different  independent 
sources  in  the  far  field,  so  the  noise  at  a  hydrophone  can  be  considered  as  the 
summation  of  a  large  number  of  uncorrelated  plane  waves  propagating  from  var- 
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ious  directions.  Such  a  field  is  usually  described  by  a  directional  density  func¬ 
tion  which  indicates  the  amount  of  noise  power  coming  from  each  direction. 
This  model  can  be  simulated  by  many  statistically  independent  point  sources 
randomly  distributed  on  the  surface  of  a  very  large  sphere  with  the  omnidirec¬ 
tional  hydrophones  located  at  the  center  of  the  sphere. 

Assume  the  noise  field  to  be  an  ergodic,  wide  sense  stationary  stochastic 
process,  and  the  noise  fields  at  different  frequencies  to  be  statistically  indepen¬ 
dent.  Then,  take  two  omnidirectional  hydrophones  separated  by  a  distance  d 
with  the  origin  located  at  the  midpoint  between  the  sensors.  Fig.  5-7.  Then,  by 
Parseval’s  theorem,  the  directional  average  noise  intensity  arriving  at  the  origin 
from  an  elemental  spherical  surface  is 


£[|«0^6-4>)|2i  =  1  |  N0(f,e,$)\2df 


and  the  total  average  noise  intensity  at  the  origin  is 


(5-8) 


£[|»o«)|2]  =  J  •  <5-9> 

47C-00 

where  E  stands  for  the  expected  value,  over  bar  stands  for  average  value,  and  d£2 
is  elemental  solid  angle. 
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FIGURE  5-7.  Polar  coordinate  system  used  in  noise  field  description. 


Let  and  x2  be  the  position  vectors  (r,0,<j>)  describing  the  location  of  the 
first  and  second  sensors,  respectively,  and  let  u  be  a  unit  vector  in  the  direction 
opposite  to  that  in  which  the  plane  wave  is  propagating.  Then  if  a  plane  wave 
propagates  past  the  sensors,  the  output  of  the  second  sensor  is  a  delayed  version 
of  the  output  of  the  first  sensor.  That  is, 

nx  (t,  0,  <}>)  =  n0(f  +  T,  6,  <f>)  (5-10) 

/z2  (r,  0,  <J))  =  nQ(t-z,  0,<J>)  (5-11) 

where 

x-^u-  (xl-x2)/c  .  (5-12) 

The  directional  cross  correlation  between  the  two  sensors  is 
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fl12(2t,e,<j>)  =  £[/i0(r(e,<i>)K0*  (f+2x,e,<j))]  .  (5-13) 

Since  the  noise  is  a  wide  sense  stationary  random  process,  then  the  noise  signal 
can  be  presented  via  Fourier-Stieltjes  integral  [25]  [26] 

oo 

n0  (t,  0,  <>)  =  j  e~j2nfldT\  (ft  0,  (?)  (5-14) 

where  dr\  (f)  is  called  the  random  amplitude.  Since  the  noise  is  a  zero  mean  ran¬ 
dom  process,  then 


(dr\(f))  =  0 

Eq.  (5-13)  can  be  written  as 


(5-15) 


e  [«0  (t,  e>  (j>)  «o*  (r + 2<r»  e,  <}>)  ]  =  J  f  (dr\(f)dT]*(f)) 

— oo— oo 

X  exp  [~j2n  (J-f)  t  +j2nf(2T)  ]  (5_16) 

Since  Eq.  (5-16)  must  be  a  function  of  t  only  (the  process  is  WSS),  then 

(dT\(f)dT]*(f))  =  S(f)b(f-f)dfdf  ;  (5-17) 

the  spectra  dr\  (f)  at  different  frequencies  must  be  uncorrelated. 

S  (/) )  is  the  spectral  density  of  the  random  noise  and  represents  the  amount 
of  “power  density”  at  different  frequencies 

S(f)  =  \N0(f,Q,<?)\2  .  (5-18) 

Therefore,  Eq.  (5-13)  can  be  written  as 


tf12(2T,e,<}>)  =  j\N0(f,e,(?)\2^2x)df 


— oo 


(5-19) 
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and  substituting  for  xfrom  Eq.  (5-12),  we  get 


*12( 2t,e,4o  =  J|N0(/,e,({))jV2n/(w  (JCl  Xl)/C)df 


Integrating  over  all  possible  directions  of  noise  arrival,  we  get 


(5-20) 


*12(2t)  =  J  J|^0(/,e,<}))|V‘2n/("‘(JCl  Xi)/c)dfdQ  .  (5-21) 

4n-oo 

Interchanging  the  order  of  integration,  we  get 


*12  (2x)  =  jrj|^o^e.4>)jV27C/(“‘ 

— oo  _4jc 


(jc  1-x2)/c) 


dCl 


df 


(5-22) 


Notice  that  R12  depends  on  the  separation  of  the  sensors  and  the  orientation 
of  the  sensors  relative  to  the  direction  from  which  the  wave  is  propagating  but 
not  the  absolute  position  of  the  sensors.  A  field  with  these  properties  is  spatially 
homogeneous.  Fields  composed  of  a  superposition  of  uncorrelated  plane  waves 
are  spatially  homogeneous:  the  cross  correlation  is  invariant  under  translation 
but  not,  in  general,  under  rotation  of  the  position  of  the  sensors. 


The  inner  integral  in  Eq.(5-22)  is  in  the  form  of  a  Fourier  transform  of  the 
noise  angular  spectral  density  at  frequency  f  to  the  linear  domain  of  separation  d 
along  the  array.  Thus,  the  single  frequency  spatial  cross-correlation  for  the  sig¬ 
nals  at  points  1  and  2  is 


Gn{d,f)  =  J|A0(/,e,<}))|V2,t/(l,'‘//c)dQ  (5-23) 

4n 

where  d  =  an<*  u  d  =  d sin<{>,  where  0  is  the  elevation  angle,  Fig.  5-7. 

Then 
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n  tt/2  ..  /sin <j) 

r  r  - o  J2nf—z- 

G12  (d,f)  =  J  J  |iV0(e,  4>,/)|^  c  cos tydtydQ  .  (5-24) 

-7C-7T/2 

Gl2(d,f)  is  called  the  spatial  cross-spectral  density  function. 

Thus,  the  time-space  correlation  function  of  the  ambient  noise  is 


Rn{d~x)  =  j  Gu(d,f)el2nfxdf  (5-25) 

— oo 

which  represents  the  Fourier  transform  of  the  spatial  cross-spectral  density  func¬ 
tion. 

5-6  HAP  in  the  Presence  of  Ambient  Noise 

The  presence  of  ambient  noise  will  affect  the  performance  of  the  HAP  algo¬ 
rithm.  It  will  dictate  a  lower  limit  on  the  strength  of  the  detectable  target  with  a 
given  signal-to-noise  ratio.  We  proceed  in  analyzing  the  effect  of  the  ambient 
noise  on  the  HAP  algorithm  using  the  following  assumptions: 

•  The  noise  is  a  spatially  homogeneous,  temporally  wide  sense  stationary  and  ergodic 
random  process; 

•  The  noise  is  a  Gaussian  random  process  with  zero  mean  and  temporal- spatial  cross-cor¬ 
relation  R  -(d,%) ,  where  i  and  j  designate  the  i-th  and  the  j-th  sensors  and  d  is  the  sep¬ 
aration  length  between  them. 

5-6.1  Target  Contaminated  with  Noise 

As  a  first  step  in  evaluating  the  effect  of  ambient  noise  on  the  performance 
of  the  holographic  array  processor  we  assume  that  only  the  target  signal  is  con¬ 
taminated  with  additive,  spatially  correlated  ambient  noise  with  mean  equal  to 
zero  and  strength  equal  to  a2  ,  Thus  the  holographic  image  at  reference  depth 
zRi  is 

h 

H  (zRi)  =  |  [Pt  (z,  rQ;zt,  rf)  +  nt  (z)  ]  PR*  (z,  rQ;zRi,  rR)  dz 
0 


(5-26) 


no 


Since  the  noise  is  an  ergodic  zero-mean  normal  process,  the  average  holo¬ 
graphic  image  is  equal  to  the  holographic  image  in  the  absence  of  the  noise 

(H(zRi))  =  H  (zRi)  (5-27) 

and  the  cross-correlation  between  the  holographic  image  at  depths  zRiand  zRj  is 

hh 

(H  (zRi)  H  (zRj)  >  =  H  (zRi)  if  (zRj)  +  JJ*n,  (U-~z\)  PRi  (z)  PR*  (i)  dzdz  (5-28) 

00 

where,  Rn{  is  the  spatial  correlation  function  of  target  noise,  and  PRi  and PRj  are 
the  pressure  fields  from  reference  sources  at  depth  zR  ■  and  zRj  as  measured  at  the 
receiving  array. 

The  reconstructed  target  field  is 

hK  * 

Q(z',r')  =  J  G(zR;z',r')H  ( zR)dzR  (5-29) 

0 

where,  G  ( zR;z r')  is  the  pressure  field  from  a  point  test  source  at  location 
(r',  z')  as  measured  (estimated)  at  the  virtual  array,  and  we  assume  that  it  is 
noise-free.  Taking  the  average  of  Eq.  (5-29),we  find  that  the  reconstructed  field 
has  a  mean  value  equal  to  that  without  noise, 

<<2(z',0>=  Q(z',r').  (5-30) 

The  second  moment  of  the  reconstructed  field  is 

hit  *  hR 

(Q  (a,-)  Q  (dj)  >  =  <  J  G  izjt&i)  H  (zR)  dzR  J  G*  (iR,dj)  H  (zR)  d~zR) 

0  o 

hf(hR  * 

=  J  jG(zj?;a|.)<ff  (zR)  H  (zR)  >G*  (zR;dj)  dzRdzR 

0  0  (5-31) 

where  ai  =  (r/,  z{) ,  the  i-th  location  of  the  scanning  test  source.  Substituting 
Eq.  (5-28)  for  the  cross-correlation  of  the  holographic  image,  the  second 
moment  of  the  reconstructed  field,  Eq.  (5-31),  becomes 


Ill 


hBh 


■R"R 


(Q  (&,•)  Q  (Zj) )  =  J  J  dzRdzRG  (zR;dt)  G*  ( zR,hj ) 
0  0 


x  (zr)H  (zr)  +  jjRnt(\z-z\)P (z;zR)P R(z;zR)  dzdz^  . 


(5-32) 


The  first  term  is  equal  to  the  reconstructed  field  in  the  absence  of  noise.  There¬ 
fore,  Eq.  (5-32)  can  be  written  as 


^R^R 


<Q(&i)Q  (&j))  =  QCaJQ 


/  * 


Cap 


J  J  °  i )  V  (zR,  zR)  G*  ( zR;dj )  dzRdzR  (5-33) 

0  0 


where 


hh 

V(ZR’  ~zr)  ~  \\Rnt^z~^)PR  (z‘>zR)pR(z,zR)dzdz  .  (5-34) 

00 

2  2 

Assuming  the  target  strength  is  0  ,  the  reference  source  strength  is  oR,  and 

2  1  * 

the  noise  strength  is  o^,  then 


and 


—  r\  ry  1/2 

<G(z',r')>  =  (ofap  Q  (z\  r') 


(5-35) 


(Q  Cad  Q  ( dj ) )  =  (ofoj)  (2  (a,)  (2*  (&/)  + 

^R^R 

Gn/°JJ  (zR#j)dzRdiR  . 

0  0 

Thus  the  variance  of  the  reconstructed  image  field  is 
o2q  =  (\Q(a)  -Q(d)\2) 

^R^R 

=  an,G*J  \G{zR,d)\{zR,~zR)(f  (zR,d)dzRdzR  , 

0  0 


(5-36) 


(5-37) 
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and  the  signal-to-noise  ratio  at  the  output  of  the  processor  is 


SNRq  = 


(a?cj)|e  (S)  |2 

\G(zR;a)v(zR,zR)(f  (zR,a)  dzRdzR 
0  0 


(5-38) 


We  notice  that  the  SNRq  is  a  function  of  the  target  to  noise  strength;  the 
magnitude  of  the  reference  source  power,  oRy  has  no  effect  on  the  SNR  at  the 
output  of  the  processor.  In  essence,  the  second  moment  of  the  reconstructed 
image  can  be  normalized  by  the  strength  of  the  reference  source  with  no  other 
consequence. 


5-6.2  Target  and  Reference  are  Contaminated  with  Noise 

The  target  and  the  reference  signals  are  contaminated  with  additive,  spa¬ 
tially  correlated  ambient  noise.  Thus  the  holographic  image  is 


h 

H  ( zRi )  =  J  {  [Pt  (z,  rQ;zr  rf;fj)  +  n(  (z;^)  ] 

0  * 

x  lPRi  r0’zRi’  rR'1 2)  +  nRi  1  >  dz  (5.39) 

where  n{  is  the  ambient  noise  associated  with  the  measurement  of  the  target  sig¬ 
nal  at  the  receiving  array,  nRi  is  the  ambient  noise  associated  with  the  signal 
from  the  reference  source  at  reference  depth  zRi  as  measured  at  the  receiving 
array,  and  tj  and  t2  are  the  times  at  which  we  measure  the  target  and  reference 
signals  respectively. 

Assume  the  two  signals  are  measured  with  a  time  interval  greater  than  the 
temporal  correlation  length  for  the  ambient  noise.  (This  correlation  length 
depends  on  the  band  width  of  the  input  filter  at  the  receiver;  see  Eq.  (5-25).) 
Therefore,  the  noise  associated  with  the  target  is  uncorrelated  with  the  noise 
associated  with  any  reference  source.  Also,  the  noise  associated  with  one  refer¬ 
ence  source  (reference  depth)  is  uncorrelated  with  the  noise  associated  with  a 
different  reference  source  (different  reference  depth),  since  the  measurements 
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from  different  reference  sources  are  taken  with  a  time  interval  greater  than  the 
temporal  correlation  length  of  the  noise.  Taking  this  into  account,  and  assuming 
that  the  noise  signal  is  zero-mean,  the  first  moment  of  the  holographic  image  is 


(H  (zRi)  )  =  J  [P  t  (z,  r0;zr  r{)  PR*  (z,  rQ;zRi,  rR)  +  P(  (z,  r0;z,,  r,)  < nRi  (z)  ) 
0 

+  PR*  ( z ,  rQ;zRi,  rR)  ( nt  (z)  >  +  <n,  (z)  nRi  (z)  >  ]  dz 


(5-40) 


where  the  average  values  in  the  integrand  are  equal  to  zero.  Therefore,  the  mean 
value  of  the  holographic  image  is 


(H(zRi))  =  H(zRi) 


Ri> 


(5-41) 


which  is  equal  to  the  holographic  image  in  the  absence  of  ambient  noise. 
The  second  moment  of  the  holographic  image  is 


hh 

(H  (zRi)  H  (zRj)  >  =  JJ  {  [Pt  (2;zt)  +  nt  (z)  ]  [ PRi  {z\zRi)  +  nRi  (z)  ]*  } 

00 


x  {  [Pt(z;zt)  +  n(  (z)  ]  [PRj(.hzRj)  +  nRj(z)]*  }*  dzdz 


(5-42) 


where  the  dependence  of  the  pressure  fields  on  the  corresponding  ranges  is 
dropped  out,  merely  for  compact  presentation.  After  multiplication,  and  some 
simplification  using  the  noise  properties,  we  get 

hh 

(H(2Ri)H  (lRj))  =  «(zR,)W*(V+JJ{(P  ,(z )P,*(i)  (nRi(z)nRjCz)> 

00 

+  PR*  (■ z\zRi )  PRj  (z;zRj)  (nt  (Z)  nt  (z) ) 

+  (nt(z)nt(z)nRi(z)nRj(z))}dzdz  . 

Note  that  ( nt  (z)  n(  (z) )  =  Rn(  (|z  -  z))  is  the  cross-correlation  function  of  the 
noise  associated  with  the  target  signal  as  measured  at  the  receiving  array.  The 
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cross-correlation  function  of  the  noise  associated  with  the  reference  sources  as 
measured  at  the  receiving  array  is 


(nRi(z)nRi(z))  =  { 


RnR  ;  *  ~  J 


■  *  o  ;  i  *j 

Since  the  ambient  noise  is  a  zero  mean  Gaussian  process,  then 

(nt(z)nt(z)nRi(z)nRj(z))  =  (n(  (z)  w,  (£)  )(nRi  (z)  nRj  (z) )  + 

(nt  (z)  (z)  ><«f  (z)  (z)  >  + 


Thus,  Eq.  (5-43)  can  be  written  as 


(5-44) 


(nt(z)  nRj(z))(nt(z)  nR-(z))  (5-45) 


(fnzRi)H  (**;)>  = 

hh 

00 

Rnt  (lz_  z\)RRi  (z’zRi)RRj(*’zRj)  + 

4  Rnt{\z-%)RnR(\z~z\)  ]dzdz  i=j 

hh 

JK  <  I2  -  Zl )  PR?  (*>ZRi)  PRj &zRj)  dzdz  i*j  . 

00  (5-46) 


Scanning  the  volume  of  the  ocean  in  the  vicinity  of  the  target  with  a  test 
point  source,  and  weighting  the  holographic  image  at  the  virtual  array  by  the 
estimated  pressure  from  the  test  source,  we  get  the  reconstructed  image  field 

kg  * 

Q(z',r')  =  \G(zR,z',r')H  ( zR)dzR  . 

0 

The  first  moment  of  the  reconstructed  field  is 


(5-47) 
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k*  _* 

(Q  (z\  r')  >  =  J  G  ( zR;z',  r')  (H  (zR)  )dzR  =  Q  ( z',  r')  (5-48) 

0 

which  is  the  same  as  in  the  previous  section.  The  mean  value  of  the  recon¬ 
structed  field  is  equal  to  the  noise-free  reconstructed  field.  The  second  moment 
of  the  reconstructed  field  image  is 

ha  kg  * 

(Q  (af)  Q  (hj) >  =  J  J  G  ( zR;ai )  < H  ( zR )  H  (zR)  )(?  Gr&P  dzRdlR  *  t5"49) 
0  0 

Substituting  for  the  second  moment  of  the  holographic  image,  Eq.  (5-46),  we  get 


< Q  (dt)  Q  (dj) >  =  Q  ( ^ )  <f  ( d} )  +  J  J  G  (z^)  v  (zR,  zR)  G*  (zR;dj)  dzRdzR( 5-50) 

0  0 


where 


hh 

00 

^(U-^)V(z;z/?)^(i;^)  + 

/?n,(lz-z|)i?rt^(lz-z|)  ]dzdz  i=; 

hh 

\\Rn,(\z-l\)PK*  (,r,zR)PR(h'zR)dzdi  .  . 

00  (5-5D 

In  comparison  with  Eq.(5-46),  we  substituted  PRi(z’,zRi)  by  (z;z/p ,  and 
PRj(z;zRj)  by  PR(z;zR) .  The  variance  of  the  reconstructed  field,  Q,  is 

o*(a)  =  (le(a)-e(a)!2)  (5-52) 

where  we  used  Eq.  (5-48)  for  the  (Q  (a) );  therefore,  Eq.  (5-52)  can  be  written  as 

Og  (a)  =  (Q  (a)  2*  (a))-fi  (a)  (2*  (a)  .  (5-53) 

Substituting  Eq.  (5-50)  in  Eq.  (5-53),  we  get 
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^R^R 

CTg(fl)  =  J  jG(zR;a)v(zR,zR)(f  (zR\a)dzRdzR  (5-54) 

0  0 

where  v(zR,zR)  is  defined  in  Eq.  (5-51). 

Given  the  temporal-spatial  correlation  function  of  the  ambient  noise,  one 
can  find  the  cross-correlation  matrix  of  the  reconstructed  field  using  Eqs.  (5-50) 
and  (5-51).  The  variance  of  the  reconstructed  field  can  be  found  by  subtracting 
the  mean  square  of  the  reconstructed  field  from  the  diagonal  of  the  cross-corre¬ 
lation  matrix,  Eq.  (5-53). 

5-6.3  Isotropic  Ambient  Noise  Model 

It  is  sufficient  to  assume  an  isotropic  model,  which  is  the  least  complex 
mathematically,  and  yet  of  practical  interest.  In  the  mid-band  frequency,  50  Hz 
to  500  Hz,  where  the  distant  shipping  and  the  surface  effects  are  the  dominating 
sources  of  noise,  the  total  ambient  noise  at  an  omnidirectional  hydrophone  has  a 
nearly  isotropic  angular  distribution. 

Fig.  5-8  (a)  shows  the  angular  intensity  of  an  isotropic  ambient  noise  model 
at  frequency  /  as  would  be  measured  by  an  omnidirectional  hydrophone,  while 
Fig.  3-8  (b)  shows  the  spatial  cross-spectral  density  function  for  this  model. 
Notice  that  when  the  separation  between  the  sensors  is  a  multiple  of  X/2,  the 
cross  spectral  density  is  equal  to  zero,  and,  as  a  consequence,  the  noise  signals 
measured  at  these  sensors  are  statistically  independent. 

2  2 

Assuming  a  constant  angular  intensity,  |  ( N0  (0,  <}>,/) )  |  =  |  (NQ  (f) )  |  ,  the 
total  spectral  noise  intensity  at  frequency  f  is 


in(f)  =  l\(.N0<f»fdn 

4n 

=  47tj  (NQ  (/)  )|2 


(5-55) 


The  cross  spectral  density,  Eq.  (5-24),  becomes 
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tc  b/2  ,„_/sin4> 

00  9 

Gl2(d,f)  =  J  J  |N0(/)p  cos<M<M9. 

-71-71/2 

Changing  variables,  Eq.  (5-56)  can  be  written  as 


7t  1/X 


c,2wy)  =  j-j  \N0(/)\y2,,v\d^ 


sin  (())) 


-7I-1/X 
c 


where  £  =  — j— -«  and  ^  ~  CarryinS  out  the  integration  we  get 


(5-56) 


(5-57) 


9sin(27td/X)  _ 1X 

an(d,f)  =  **p0(f>\2  ■ .  (5-58) 

therefore 


Gn(d,f)  =  47C|iV0 (/) |2sinc (2rf/X.)  .  (5-59) 


FIGURE  5-8.  (a)  Angular  intensity  distribution  of  the  isotropic  ambient  noise;  (b) 
Normalized  cross- spectral  density  for  the  isotropic  ambient  noise. 
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Thus  the  time-space  cross  correlation  function  for  band  limited  noise  is 


Rn(d,x)  =  j  4K\N0(f)\2sinc(2d/'k)W(f)ei2nfxdf 


(5-60) 


where  W(f)  is  a  weighting  window  characterizing  the  frequency  band  of  interest. 
Assuming  that  the  spectral  density  is  constant  over  the  band  of  interest, 

|Nq  (J)  |  =  A ,  and  the  weighting  window  is 

1  for  /  •  < /<  fnr 

■  *  J  min  J  J  max 


W(f)  =  { 


0 


else 


then 


(5-61) 


/- 


Rn(d,x)  =  J  47iAsinc  (2d/\)  e'2n^df  . 

f wi* 

For  a  time  delay  T  =  0,  the  spatial  correlation  function  is 


(5-62) 


Rn(d,0)  =  4  nA  J  sine  (2d/'k)  df 


L 


471/4 


Sl( — )  -Si  ( — - — ) 


2d  L 


where  the  special  function,  sin  integral  [27],  is 


(5-63) 


( sin  (x) 

•Si  (z)  =  J — - — dx 


(5-64) 


For  more  detail  about  this  function  see  appendix  II. 


5-6.4  Simulation  Results 

Figs.  5-9,  5-10,  and  5-11  show  the  magnitude  of  the  normalized  mean  of  Q, 
solid  lines,  and  the  effect  of  adding  one  standard  deviation  to  the  mean  value, 
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dotted  lines,  for  various  normalized  (transmission  loss  is  not  included)  target 
strength.  The  ambient  noise  is  considered  to  be  isotropic  with  average  intensity 
of  70  dB.  Fig.  5-9  presents  these  plots  for  a  normalized  reference  source  power 
of  100  dB,  and  the  normalized  target  power  varies  from  90  dB  to  160  dB.  Fig.  5- 
10  shows  the  same  plots  for  a  normalized  reference  source  power  of  110  dB,  and 
the  normalized  target  power  varies  from  90  dB  to  160  dB.  Fig.  5-11  shows  the 
same  plots  for  a  normalized  reference  source  power  of  120  dB,  and  the  normal¬ 
ized  target  power  varies  from  90  dB  to  160  dB.  In  all  cases  we  notice  that 
increasing  the  target  power  beyond  a  certain  point,  close  to  the  reference  source 
power,  will  not  improve  the  performance  (the  focus  and  ambiguity)  of  the  algo¬ 
rithm.  The  standard  deviation  saturates;  it  becomes  a  constant  as  a  function  of 
target  strength.  When  the  target  power  is  below  the  reference  source  power,  the 
algorithm  is  very  dependent  on  the  target  strength;  the  standard  deviation  gets 
smaller,  the  focus  gets  sharper,  and  the  ambiguity  (side  lobes)  gets  weaker  with 
increasing  target  strength. 

A  convenient  measure  of  system  performance  is  the  signal-to-noise  power 
ratio 


(5-65) 


where  Sq  is  the  noise-free  reconstructed  image  signal,  and  Nq  is  noise  associ¬ 
ated  with  the  reconstructed  image  signal. 


<|Sq|2>  =  Q{ f  (5-66) 

<|NG|2>  =  1  \G{zR’,ai)v(<zR,zR)(f  ( zR;bj)dzRdzR  . 

0  0 


(5-67) 


(C)  (0 

FIGURE  5-9.  The  mean  magnitude  of  the  reconstructed  field  (solid  line)  and  the  mean 
magnitude  plus  one  standard  deviation  of  the  reconstructed  field  (dotted  line).  The 
normalized  reference  strength  is  100  dB,  the  ambient  noise  intensity  is  70  dB,  and  the 
normalized  target  strength  varies:  (a)  90  dB;  (b)  95  dB;  (c)  100  dB;  (d)  110  dB;  (e)  120 
dB;  (f)  160  dB. _ 
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FIGURE  5-10.  The  mean  magnitude  of  the  reconstructed  field  (solid  line)  and  the  mean 
magnitude  plus  one  standard  deviation  of  the  reconstructed  field  (dotted  line).  The 
normalized  reference  strength  is  110  dB,  the  ambient  noise  intensity  is  70  dB,  and  the 
normalized  target  strength  varies:  (a)  90  dB;  (b)  100  dB;  (c)  110  dB;  (d)  120  dB;  (e)  140 
dB;  (f)  160  dB.  _ 
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(C)  (f) 

FIGURE  5-11.  The  mean  magnitude  of  the  reconstructed  field  (solid  line)  and  the  mean 
magnitude  plus  one  standard  deviation  of  the  reconstructed  field  (dotted  line).  The 
normalized  reference  strength  is  120  dB,  the  ambient  noise  intensity  is  70  dB,  and  the 
normalized  target  strength  varies:  (a)  90  dB;  (b)  100  dB;  (c)  110  dB;  (d)  120  dB;  (e)  140 
dB;  (f)  160  dB.  _ 
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Acceptable  system  performance  is  possible  if  and  only  if  the  signal-to-noise 
mean  power  ratio  is  large,  say  SNR  >  5.  Thus  the  presence  of  noise,  for  a  given 
reference  source  power  level,  will  impose  a  lower  limit  on  the  detectable  target 
signal.  From  the  above  analysis,  we  notice  that  the  noise  part  of  the  recon¬ 
structed  field  is  a  function  of  both  the  target  and  reference  signals.  Therefore, 
for  a  given  reference  source,  increasing  the  target  power  level  will  lead  to 
increasing  both  the  signal  and  the  noise  terms  of  the  reconstructed  field;  this  in 
turn  leads  to  the  conclusion  that  increasing  the  target  power  beyond  a  certain 
level  will  lead  to  no  improvement  in  the  SNR. 

Fig.  5-12  shows  the  SNR  at  the  output  as  a  function  of  target  strength  nor¬ 
malized  by  the  transmission  loss,  for  various  normalized  reference  source 
strengths,  as  shown  on  the  curves.  We  notice  that  the  SNR  is  a  weak  function  of 
target  strength  for  target  power  greater  than  the  reference  source  power,  and  it  is 
a  strong  function  of  the  target  strength  when  the  latter  is  less  than  the  reference 
source  power. 

Transmission  loss  due  to  wave  propagation  in  the  ocean  was  reviewed  in 
chapter  1.  As  an  example,  assume  the  goal  is  to  detect  a  target  at  range  100  km 
with  minimum  SNR  of  15  dB,  given  the  reference  source  power  is  190  dB.  Since 
the  reference  source  is  near  the  target,  then  the  signals  from  both  will  endure 
approximately  the  same  amount  of  transmission  loss.  Assume  the  transmission 
loss  is  80  dB;  then  the  normalized  reference  source  power  is  110  dB.  From  the 
plot  and  for  SNR  equal  to  15  dB,  we  find  the  normalized  target  strength  to  be 
102  dB,  so  the  minimum  detectable  source  strength  (with  the  given  minimum 
SNR)  is  182  dB.  Therefore,  a  target  at  that  range  (transmission  loss)  with  power 
of  182  dB  or  greater  will  be  detected  with  SNR  of  15  dB  or  better. 


124 


FIGURE  5-12.  Signal  to  noise  ratio  of  the  reconstructed  field  as  a  function  of  normalized 
target  strength,  and  for  various  normalized  reference  source  power  levels  (as  shown  on  the 
curves).  The  ambient  noise  intensity  is  70  dB. _ 


GENERAL  CONCLUSIONS 


This  dissertation  has  demonstrated  through  theoretical  analysis  and  numeri¬ 
cal  modeling  that  the  holographic  array  processor  can  focus  with  great  accuracy 
the  acoustic  field  from  a  point  source,  in  a  range  and  depth  dependent  environ¬ 
ment. 

It  has  been  shown  that  truncating  the  virtual  or  the  receiving  array  reduces 
the  imaging  performance  of  the  holographic  array  processor.  Analytical  and 
numerical  analyses  showed  that  the  processor  performance  is  more  severely 
degraded  by  truncating  the  virtual  array. 

The  performance  of  the  holographic  array  processor  in  the  presence  of  ambi¬ 
ent  noise  has  been  analyzed.  Theoretical  analysis  and  numerical  simulation 
revealed  that  the  performance  is  limited  by  the  signal  to  noise  ratio  of  the  refer¬ 
ence  source  signals. 

The  following  areas  are  suggested  for  further  research: 

•Multiple  target  detection  and  tracking  using  the  HAP  method. 

•The  effect  of  a  mismatch  in  reference  and  target  frequencies. 

•The  performance  of  the  HAP  in  the  presence  of  internal  waves. 

•The  performance  of  the  HAP  under  nonadiabatic  propagation  conditions. 

•The  performance  of  the  HAP  in  shallow  water. 
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Appendix  I 


WKB  Approximation 

The  WKB  approximation  is  used  to  calculate  an  approximation  to  the  eigen- 
numbers.  Because  the  WKB  approximation  is  valid  in  the  geometric-acoustic 
limit,  the  mode  expansion  presented  here  is  simply  related  in  many  ways  to  a  ray 
theory  solution.  The  depth  dependent  portions  of  the  field  are  standing  waves 
that  must  satisfy  the  differential  equation 
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where  the  vertical  propagation  constant  is  defined  by 

im  =  k2(z’r)-^m- 
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FIGURE  1-1.  Outward  radiating  cylindrical  wave _ 

The  mode-propagation  wavenumbers  are  determined  from  the  trapping 
requirement  that  the  traveling  waves  that  constitute  each  mode  must  be  out  of 
phase  by  m7l  radiance. 

The  WKB  approximation,  or  phase  integral  method,  yields  solutions  to  Eq. 
(I- 1)  that  are  valid  far  from  the  turning  point,  or  if 
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The  WKB  approximation  yields  upward  and  downward  traveling  waves 
whose  phase  is  plus  or  minus 

2 

e(z,z\5)  =  jY(U)<iC-  (1-4) 

2* 

As  the  waves  are  out  of  phase  by  tc/2  rad  upon  reflection,  the  resonance  condi¬ 
tion  is  given  by 

Gm(zu,  z/)  =  (m  +  ^)n  (1-5) 


where  zu  and  z\  are  the  upper  and  lower  turning  points.  Referencing  the  field  to 
the  upper  turning  point,  the  depth  dependent  field  is  approximated  by 
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For  z  below  the  turning  point,  the  approximate  field  is  found  by  substituting  zl 
for  zu  and  interchanging  z  and  zl  in  the  damping  term.  This  solution  is  a  good 
approximation  when  k\z-zu\»\  and  when  the  turning  points  are  widely 
spaced.  For  z  meeting  these  conditions  such  that  0  »  1 ,  we  can  further  approxi¬ 
mate  <J>=0  for  z  above  or  below  the  turning  points,  and  as  a  result  the  ray  exit 
angle  <}>  is  very  small.  The  mode  propagation  wavenumber  can  be  written  as 

£„<r>  =  ^Tcos[<p«(r)1'  a'7> 


Using  this  notation  and  that 
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we  can  write 
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Y2  (z,  r)  = 
*m  v 


0)‘ 


r  W2(z,  r) 
N2(z„,r) 


-  cos2  cp  (r) 


c  2  (zfl,  r) 


and 


Z  r 

(0  f 

Vz-z«>  = 


<£&  -cos2^(r) 


1/2 


d-9) 


(MO) 


Appendix  II 
Sine  Integral 


The  sine  integral  is  a  special  function  and  is  defined  as 
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In  series  expansion  the  sine  integral  is 
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and  in  terms  of  auxiliary  functions  is 


7t 


Si  (z  ) )  =  ^-/(z)cosz-g(z)  sinz 


01-3) 


where  /(z)  and  g  (z)  have  asymptotic  expansion  as 
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A  rational  approximation  follows,  for  1  <  z  <  °° 
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g  (z)  = 


1  (Z4  +  c1z2  +  C2V 

z2  vz4  +  d1z2  +  d2J 


e,(z) 


(n-7) 


where  ax  =  7.241163,  a2  =  2.463936,  ^  =  9.068580,  fc2  -  7.157433 

,  =  7.547478,  c2  =  1.564072,  dj  =  15.723606,  d2  =  12.723684  . 

j  z  ^  —4 

This  approximation  has  the  error  terms  |Ej  (z)|  <  2x10  ,  and|e2(z)|  <  10 
A  more  accurate  approximation  is 


/(z) 


1 

z 


rz>  +  fljZ6  +  a2z  +  fljZ2  + 

~  ~  4  I 

^  +  bxz  +b2z  +  b3z  +l>4; 
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g(z) 

zz 


/  z8  +  CjZ6  +  C2z4  +  C3z2  +  C4  ^ 
4 


z  +  djZ  +  d2z  +^3Z  +^4 


+  e2(z) 


01-9) 


where  a,  =  38.027264,  a2  =  265.187033,  «3  =  335.677320,  a4  =  38.102495 

fe,  =  40.021433,  i>2  =  322.624911,  i>3  =  570.236280,  bt  =  157.105423 

c,  =  42.242855,  c2  =  302.757865,  c3  =  352.018498,  c4  =  21.821899 

a  =48  196927,  d2  =  482.485984,  d 3  =  1114.978885,  d4  =  449.690326 

1  _7 

This  approximation  has  the  error  terms  |Ej  (z)|  <  5x10  7,  and |e2  (z)  |  <  3x10 
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Appendix  III 

Parabolic  Equation  (PE)  Approximation  to  the  Wave  Equation 

The  propagation  of  monochromatic  sound  in  the  ocean  is  described  by  the 
Helmholtz  equation 


V2P  + 


(HI-1) 


In  cylindrical  coordinates  this  becomes 


lJL(r— ) 


2  2 

+  4-—^  +  +  fa1  (r,  z,  cp)  P 

rW  dz2 


0 


(HI-2) 


where  yi  =  - -  is  the  noroicilizcd  index  of  refraction, 

c  (r,  z,  cp)  .  „  ,  , 

ence  wavenumber,  and  c q  is  the  reference  sound  speed. 


kc  =  —  is  the  refer- 
0  c0 


Rays  propagating  with  angles  greater  than  the  maximum  exit  angle  (the 
angle  the  ray  makes  with  the  horizontal)  are  not  trapped  in  the  ocean  and  thus  do 
not  contribute  to  the  far  field.  Therefore,  for  waves  propagating  in  the  direction 
close  to  the  horizontal,  a  solution  to  the  wave  equation  may  be  attempted  in  the 
form 


P  (r,  z,  cp)  =  U  (r,  z,  cp)  H^l)  ( kQr )  (m“3) 

where  H^l)  is  the  Hankel  function  of  the  first  kind  of  order  zero  (it  describes 
the  part  of  the  solution  rapidly  varying  in  the  r  direction),  and  U  is  a  slowly 
varying  function  in  all  three  directions.  For  kQr  »  1 ,  the  Hankel  function  has  the 
asymptotic  presentation 

2  1/2 

"c”  (V)  -  exp  [i*°rl  ■ 


(ni-4) 


(ni-5) 


Substituting  Eq.  (III-3)  in  Eq.  (III-2)  and  considering  Eq.  (HI-4)  we  get 

^  +  2 ik0f  +  +  ^  ^  (n2  -  1 )  U  =  0. 

3  2  °3r  r23  2  JZ2 


The  variation  of  U  over  a  distance  of  the  order  of  the  wavelength  is  small, 
3««  «  Therefore,  the  first  term  in  Eq.  (III-5)  is  very  small,  at  least  compared 

dr  0  , 

to  the  second  term;  then 


2ik0 57  + 


i  a  2u 

r2dy2 


■  ^  +  k20(n2-l)U  =  0. 
dz2 


(HI-6) 


At  great  ranges  from  the  source  we  may  neglect  the  curvature  of  the  wavefront 
and  set 


rd<p  =  dy. 


(m-7) 


Then  Eq.  (III-6)  becomes 


,^  +  !^  +  i|  +  4(n2-  1)1/ -0. 

dy2  dz 


(HI-8) 


If  the  variation 


of  the  sound  speed  with  azimuth  can  be  neglected,  then 


2  ik 


(ni-9) 


dz 


The  pressure  at  a  point  in  the  far  field  can  be  found  by  evaluating  the  field  at 
an  initial  range1  rQ  by  some  other  method  (geometric  acoustic,  normal 
mode, ...etc.);  then  numerically,  step-by-step,  we  find  the  field  at  the  desire 

point  using  Eq.  (III-9). 


1.  r„  must  be  small  enough  that  the  field  can  be  evaluated  by  other  methods,  but  sufficiency  large  that  the 
parabolic  approximation  is  valid. 
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A  split  step  algorithm  [28]  computes  the  complex  pressure  as  a  function  of 
depth  at  a  certain  range  given  the  pressure  as  a  function  of  depth  at  a  smaller 

range; 


U(r  +  Ar,z)  =  F~l  { exp 


-iArs 


.21 


2k 


0  J 


F  [U  (r,  z)  ]  }  exp 


9  A  r 

ikQ(n2-  1)  -y 


(HI-10) 


where  r  is  range,  Ar  is  PE  range  step,  z  is  depth,  =  —  is  reference  wave- 
number,  f  is  frequency,  c0  is  reference  sound  speed,  and  5  =  kQ sincp.  The  Fourier 
transform  notation  is  used  to  denote  a  discrete  sine  transformation  from  the 
physical  space  (z-space)  to  the  vertical  wavenumber  space  (k-space). 

The  transform  is  most  efficiently  computed  for  a  number  of  mesh  points 
equal  to  integer  power  of  two.  The  set  of  depth  mesh  on  which  the  field  will  be 
defined  is  chosen  to  be  2^  equally  spaced  depths.  These  points  extend  from  the 
surface  of  the  ocean  to  well  below  the  water-sediment  interface; 


Zmax  =  2^AZ- 
max 


(HI-11) 


Once,  the  depth  mesh  is  defined,  the  k-space  mesh  is  determined; 

As  =  kQ  Asintp  (HI-12) 


where 


Asintp 


smcp 


max 


(HI- 13) 


Appendix  IV 

Wide  Angle  PE  Code  (FEPE) 

The  parabolic  equation  (PE)  method  is  an  efficient  and  accurate  approach 
for  solving  underwater  acoustic  propagation  problems  in  range-dependent  envi¬ 
ronments.  The  PE  method  replaces  the  reduced  wave  equation,  which  is  a 
boundary  value  problem,  with  an  initial  value  problem.  The  narrow  angle  PE, 
which  is  accurate  for  propagation  angles  up  to  15  degrees,  gives  rise  to  phase 
errors  that  increase  with  range  [13].  The  phase  error  problem  was  reduced  with 
the  discovery  of  the  wide-angle  PE  [29],  which  is  accurate  for  propagation 
angles  up  to  about  40  degrees.  The  phase  error  problem  was  eliminated  entirely 
with  the  discovery  of  a  family  of  higher-order  PE’s  [31]  [32],  which  are  accurate 
for  propagation  angles  of  nearly  90  degrees. 

The  acoustic  pressure  P  is  assumed  to  satisfy  the  pressure  release  boundary 
condition  P=0  at  the  ocean  surface,  the  outgoing  radiation  condition  at  infinity, 
and  the  reduced  wave  equation 

pV*  (ivP)  +K2P  =  — 4ji8 (x-xs)  CIV-1) 

where  the  complex  wavenumber  K  =  k  +  iG$\k\  is  used  to  account  for  sediment 
loss.  Absolute  value  is  used  so  that  energy  loss  occurs  in  the  direction  of  propa¬ 
gation.  The  wavenumber  is  k  —  to /c,  O  =  [407tloge]  ,  P  is  the  attenuation  in 
decibels  per  wavelength  (dBA),  p  is  the  density,  c  is  the  sound  speed,  and  x,  is 
the  source  location. 

In  cylindrical  coordinates  Eq.  (IV- 1)  becomes  [32] 

=  ik0Jl+XW  (IV-2) 


where 
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-k2  +  ^~ 

"a*2 


i5pa_^ 

p dzdz J 


OV-3) 


where  we  define  W  =  r*/2P  if  the  environment’s  range  dependence  is  weak.  To 
reduce  the  propagation  problem  to  two  spatial  dimensions,  the  azimuthal  varia¬ 
tion  was  assumed  to  be  negligible.  Also,  variation  in  range  was  assumed  to  be 

dp  ,  . 

sufficiently  weak  that  ^  can  be  ignored. 

Parabolic  equations  are  based  on  approximations  of  the  square  root  function 
appearing  in  Eq.  (IV-2).  The  narrow  angle  PE  is  derived  from  the  first  two  terms 
of  the  Taylor  series  expansion.  The  wide-angle  PE  is  derived  from  a  rational 
approximation.  Bamberger  [30]  used  Pade  series  to  approximate  the  square  root 
in  Eq.  (IV-2)  and  derived  the  following  higher-order  PE: 


Wn 

dr 


=  ik, 


V  ./> n 
^  1  +  b,  n. 
j=  1  )'n 


xu »• 
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where 


ai>  n  2n  +  1  Sm 
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L2  n  + 
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b-  =  cos 
J>n 


rJL.1 

L2/I+IJ  ’ 


and 


W  =  Wn  =  Unexp  ( ikQr )  ■ 


(IV-5) 


dV-6) 


(IV-7) 


Eq.  (IV-4)  is  referred  to  as  PEn. 

The  computer  code  FEPE  solves  PEn  using  finite  elements  for  depth  discret¬ 
ization  and  Crank-Nicolson  integration  in  range  [33].  The  tridiagonal  system 
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solver  in  FEPE  has  been  designed  to  minimize  computation  time.  The  code  uses 
an  elimination  scheme  that  involves  sweeping  downward  to  the  row  correspond¬ 
ing  to  the  ocean  bottom  to  eliminate  entries  below  the  main  diagonal  and  sweep¬ 
ing  upward  to  the  ocean  bottom  to  eliminate  entries  above  the  main  diagonal 
followed  by  back  substitution  sweeping  up  and  down  from  the  ocean  bottom.  In 
contrast  Gaussian  elimination  involves  sweeping  downward  to  eliminate  all 
entries  below  the  main  diagonal  followed  by  back  substitution  sweeping  upward. 

For  problems  involving  range-dependent  ocean  depth,  the  new  scheme  is 
more  efficient  than  Gaussian  elimination.  In  the  decomposition  into  upper  and 
lower  triangular  matrices  of  Gaussian  elimination,  it  is  necessary  to  repeat 
sweeping  downward  from  the  ocean  bottom  as  the  ocean  depth  varies.  With  the 
new  scheme,  it  is  necessary  to  repeat  sweeping  only  for  a  few  rows  near  the 
ocean  bottom.  Since  multiplication  is  faster  than  division  on  computers,  the  trid¬ 
iagonal  system  solver  has  also  been  improved  by  replacing  divisors  with  factors. 
FEPE  comes  with  a  set  of  five  starters  to  calculate  the  initial  field.  Usually  the 
Gaussian  PE  starter  or  Green’s  wide-angle  PE  starter  is  used  to  approximate 
W(0,z). 
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